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THE  NOGAPS  TEN  YEAR  AMIP  INTEGRATION 


1.  INTRODUCTION 

Numerical  global  models  of  the  atmosphere  have  reached  a  level  of  sophistication  that 
makes  them  our  most  powerful  research  tool  for  studying  large-scale  atmospheric  processes. 
These  general  circulation  models  (GCMs)  are  widely  used  for  numerical  weather  prediction 
(NWP)  and  for  climate  simulation.  Until  recently  NWP  models  had  relatively  crude  physical 
parameterizations  compared  to  climate  models,  but  the  almost  universal  adoption  of  global  NWP 
models  and  increased  emphasis  on  medium  range  prediction  (5-10  days)  has  motivated  improve¬ 
ment  of  NWP  model  physics.  Now  there  is  usually  little  difference  between  GCMs  used  for 
NWP  and  for  climate  studies  except  that  operational  NWP  models  are  normally  run  with 
significantly  higher  resolution  than  are  climate  GCMs. 

There  are  now  dozens  of  GCMs  being  run  at  atmospheric  research  organizations  and 
operational  NWP  centers  around  the  world.  The  complexity  of  atmospheric  processes  has  lead 
to  a  considerable  lack  of  unanimity  among  atmospheric  scientists  about  how  to  model  these 
processes.  Therefore  there  is  great  diversity  among  the  GCMs'  physical  parameterizations.  Such 
diversity  is  valuable  as  a  path  toward  determining  the  "best"  parameterizations,  provided  there 
is  some  kind  of  rigorous  model  intercomparison  effort  under  controlled  conditions. 

The  U.S  Department  of  Energy's  Lawrence  Livermore  National  Laboratory  (LLNL)  has 
recognized  the  need  for  comprehensive  GCM  intercomparison  and  established  the  Program  for 
Climate  Model  Diagnosis  and  Intercomparison  (PCMDI).  In  cooperation  with  the  World  Climate 
Research  Programme's  Working  Group  on  Numerical  Experimentation  (WGNE)  PCMDI  has 
established  the  Atmospheric  Model  Intercomparison  Project  (AMIP)  as  a  way  of  bringing  many 
of  the  world's  GCMs  together  for  a  comprehensive  intercomparison.  The  numerical  experiment 
run  by  all  participants  is  a  10  year  integration  over  the  years  1979-1988,  with  PCMDI  providing 
prescribed  sea  surface  temperature  (SST)  and  sea  ice  coverage  for  this  period.  The  solar  constant 
and  atmospheric  C02  concentration  is  also  specified  by  PCMDI.  All  other  model  parameters  and 
degrees  of  freedom  such  as  resolution  are  at  the  discretion  of  the  participants. 

Participation  by  the  world's  atmospheric  modeling  groups  in  AMIP  has  been  almost 
unanimous.  AMIP  is  by  far  the  largest  numerical  experiment  of  this  kind  ever  attempted.  In 
addition  to  basic  research  modeling  groups,  every  major  operational  NWP  center  in  the  world  is 
participating,  clear  evidence  that  these  centers  believe  their  models  are  legitimate  climate  GCMs 
and  can  be  evaluated  as  such.  PCMDI  provides  data  management  and  other  logistical  support 
to  the  participating  GCM  groups,  and  is  also  supporting  several  groups  with  computer  time  on 
LLNL  computer  systems.  PCMDI's  ultimate  goal  is  to  collect  the  10  year  model  output  histories 
from  each  participating  group.  They  are  requesting  each  group  to  output  histories  every  six  hours 
of  integration,  so  the  volume  of  data  to  be  handled  is  enormous;  about  2  terabytes. 
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The  Navy  Research  Laboratory  (NRL)  is  participating  in  AMIP  because  the  Navy  has 
great  interest  in  improved  global  atmospheric  forecasts.  The  Navy  Operational  Global 
Atmospheric  Prediction  System  (NOGAPS),  developed  over  the  past  several  years  by  Navy 
atmospheric  modelers  in  Monterey,  Ca.,  generates  or  directly  influences  a  wide  variety  of 
meteorological  and  oceanographic  operational  products  for  Navy  fleet  users.  NOGAPS 
improvements  will  directly  benefit  these  users  by  improving  every  NOGAPS  product.  AMIP  is 
an  opportunity  to  compare  NOGAPS  to  the  other  major  GCMs  around  the  world,  showing 
strengths  and  weaknesses  that  can  guide  research  on  future  model  improvements.  The  10  year 
AMIP  NOGAPS  history  is  also  an  invaluable  data  set  supporting  a  wide  variety  of  NRL  research 
projects. 

A  summary  of  the  NOGAPS  forecast  model  features  and  some  details  of  how  the  model 
was  modified  for  the  AMIP  experiment  are  given  in  section  2  of  this  report.  The  large  volume 
of  model  history  data  produced  by  the  NOGAPS  AMIP  integration  required  special  attention  to 
data  organization  and  storage.  Section  3  describes  the  data  management  strategy  used  for  the 
NOGAPS  AMIP  data  set.  In  section  4  the  model  parameters  available  in  the  AMIP  data  set  are 
itemized.  Concluding  comments  and  data  access  procedures  are  given  in  section  5.  The 
appendix  contains  listings  of  Fortran  software  which  is  available  to  facilitate  access  to  the  AMIP 
data  and  post-processing  of  the  raw  history  files. 

2.  NOGAPS  AMIP  MODEL  DESCRIPTION 

A  detailed  description  of  the  NOGAPS  forecast  model  is  in  Hogan  et  al.  (1992).  Table  1 
summarizes  the  major  features  of  the  model.  Briefly,  the  model  is  a  global  spectral  model  similar 
in  concept  to  the  operational  global  NWP  models  run  by  other  large  operational  centers.  The 
physical  parameterizations  are  representative  of  those  used  in  other  NWP  models  and  climate 
GCMs  run  by  many  major  meteorological  research  centers.  The  version  of  NOGAPS  run  for 
the  AMIP  10  year  integration  is  very  similar  to  the  system  run  operationally  by  Fleet  Numerical 
Oceanography  Center  (FNOC)  in  Monterey.  There  are  no  significant  differences  be'ween  the 
physical  parameterizations  of  the  operational  NOGAPS  and  the  AMIP  model.  The  AMIP 
integration  was  done  on  the  CRAY  YMP  at  Stennis  Space  Center  (SSC),  Micsissippi,  so  the 
floating  point  results  are  64  bit  precision,  while  the  operational  NOGAPS  runs  on  the  FNOC 
CYBER  205  using  32  bit  operations.  No  meteorologically  significant  differences  have  been 
observed  because  of  this  word  size  difference,  however.  To  keep  computer  time  costs  within 
reasonable  bounds  the  AMIP  NOGAPS  has  a  horizontal  resolution  of  triangular  truncation 
wavenumber  47  (T47)  rather  than  the  T79  resolution  of  the  operational  NOGAPS.  The  number 
of  levels  is  18,  the  same  as  used  operationally.  Interestingly,  T47  is  a  significantly  higher 
resolution  than  most  of  the  other  participating  AMIP  models. 
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Table  1:  Description  of  NOGAPS  Forecast  Model. 


-  SPECTRAL  TRIANGULAR  TRUNCATION 

-  RESOLUTION 

CURRENT  OPERATIONAL:  T79/L18 
FUTURE:  -T150/L36 

COMPUTATIONAL  DETAILS 

-  SEMI-IMPLICIT  TIME  DIFFERENCING 

-  IMPLICIT  ZONAL  ADVECTION 

VORTICITY 

MOISTURE 

-  SILHOUETTE  OROGRAPHY 

-  PLANETARY  BOUNDARY  LAYER 

EXPLICITLY  RESOLVED 

STABILITY  DEPENDENT  MIXING  COEFF. 

SHALLOW  CUMULUS 

PREDICTED  GROUND  HYDROLOGY 

-  RADIATION 

DIURNAL  CYCLE 

PREDICTED  STABLE/CUMULUS  CLOUDS 
OZONE  HEATING 

-  CUMULUS 

ARAKAWA-SCHUBERT 

DOWNDRAFTS 

MASS  FLUX  KERNEL  OR  "RELAXED"  SCHEMES 

-  LARGE  SCALE  PRECIPITATION 


-  GRAVITY  WAVE  DRAG 


A  mean  solar  constant  of  1365  watts  m'2  and  a  global  mean  C02  concentration  of 
345  ppm  was  specified  by  PCMDI  for  all  participating  models.  These  are  standard  input 
parameters  of  the  model  and  are  quite  close  to  the  values  used  in  the  operational  NOGAPS.  The 
bottom  boundary  conditions  of  SST  and  sea  ice  required  some  model  modification,  however. 
When  NOGAPS  runs  as  an  operational  NWP  model,  the  SST  and  sea  ice  is  updated  at  the 
beginning  of  each  forecast,  but  held  fixed  during  typical  NWP  integrations  of  a  week  or  less. 
For  the  ten  year  AMIP  integration  NOGAPS  was  modified  to  allow  continuous  updating  of  the 
SST  and  sea  ice,  as  is  done  in  other  climate  GCMs. 

As  stated  above,  PCMDI  provided  the  SST  and  sea  ice  coverage  to  be  used  by  all  AMIP 
models.  These  are  monthly  mean  2.0  x  2.0  deg  fields  covering  the  1979-1988  period  and  were 
prepared  jointly  by  Climate  Analysis  Center  (CAC)  of  NOAA  and  the  Center  for  Ocean-Land- 
Atmosphere  (COLA)  Interactions  at  the  University  of  Maryland.  The  monthly  averages  are 
assumed  valid  at  the  middle  of  each  month.  The  SST  and  sea  ice  fields  are  bi-cubic  spline 
interpolated  to  the  approximately  2.5  deg  Gaussian  transform  grid  of  the  T47  NOGAPS  and 
linearly  interpolated  in  time  every  model  time  step  to  ensure  smoothly  varying  bottom  boundary 
conditions. 

3.  AMIP  DATA  MANAGEMENT  STRATEGY 

A  numerical  experiment  on  the  scale  of  the  10  year  NOGAPS  AMIP  integration  requires 
careful  consideration  of  the  organization  of  the  model  history  files  to  allow  reasonably  flexible 
access  to  the  data.  The  NOGAPS  AMIP  history  is  about  100  Gbytes,  so  it  presents  a  daunting 
challenge  to  anyone  interested  in  exploiting  this  unique  data  set.  Special  software  was  written 
for  the  forecast  model  to  allow  easy  model  restarts  and  reruns  as  needed  during  the  integration. 
This  software  is  also  useful  for  model  history  post-processing.  The  Fortran  subroutines  are  listed 
in  the  appendix  and  are  available  on  the  SSC  computer  system. 

To  understand  the  data  management  philosophy  used  for  the  NOGAPS  AMIP  model 
integration  history,  a  brief  description  of  the  SSC  computer  system  is  necessary.  The  heart  of 
the  system  is  an  8  processor,  128  Mword  CRAY  YMP  with  about  40  Gbytes  of  on-line  disk 
storage.  A  second  single  processor  YMP  acts  as  a  dedicated  file  server  connected  to  the  large 
machine  via  a  1.0  Gbit/sec  network.  The  file  server  also  has  about  40  Gbytes  of  disk  storage  and 
is  connected  to  a  tape  cartridge  archive  system  with  several  Terabytes  of  capacity.  Several  other 
smaller  computers  and  workstations  are  also  on  the  high  speed  network.  The  important 
distinction  between  the  two  CRAY'S  is  the  way  the  disk  storage  is  managed  on  the  two  systems. 
The  big  machine's  disks  are  all  defined  as  "scratch"  space,  with  file  duration  of  only  a  few  hours 
or  days,  depending  on  file  size  and  access  frequency.  Access  to  these  scratch  files  is  direct, 
however,  so  I/O  performance  is  quite  good. 
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The  disks  physically  attached  to  the  file  server  CRAY,  on  the  other  hand,  are  defined  as 
a  "shared"  file  system  among  all  the  machines  on  the  high  speed  network,  including  the  big 
CRAY.  This  means  that  jobs  running  on  the  big  CRAY  can  functionally  access  files  on  this 
shared  file  system  just  like  scratch  files.  There  is  a  big  performance  penalty,  however,  since 
shared  file  data  must  be  physically  moved  from  the  file  server  to  the  big  machine  across  the  high¬ 
speed  network.  A  high-speed  network,  even  with  a  1.0  Gbit  bandwidth,  is  much  slower  than 
direct  access  disks. 

Files  on  the  shared  system  are  subject  to  migration  to  the  archive  system  under  the  control 
of  sophisticated  UN1COS  file  management  software,  also  based  on  file  size  and  access  frequency. 
Retrieval  of  files  from  the  archive  is  "transparent"  to  users,  except  for  noticeable  time  delays  as 
the  data  moves  from  the  archive  tape  cartridges  back  to  the  shared  disks.  During  heavy  system 
use  periods  these  delays  can  be  quite  significant  as  network  traffic  overwhelms  the  ability  of  the 
system  to  efficiently  move  data. 

The  design  of  this  multi-level  file  system  structure  dictated  that  some  care  be  taken  in 
running  data  intensive  computations  such  as  the  NOGAPS  AMIP  integration.  The  model 
typically  ran  one  month  integration  period,  during  which  model  history  was  written  to  "scratch" 
files.  At  the  end  of  each  month's  integration,  this  history  was  copied  to  a  shared  file  system 
directory  using  the  UNIX  "tar"  utility.  In  the  process  about  500  history  files  were  consolidated 
into  5  tar  files  and  2  small  auxiliary  files. 

Table  2.  Example  of  AMIP  directory  "pcmdi0885",  showing  5  tar  files  and  2  small  auxiliary 
files.  The  expanded  contents  of  the  tar  files  are  shown  in  Table  3. 


-rw-rw-r—  1  rosmfnoc 
-rw-rw-r—  1  rosmfnoc 
-rw-rw-r--  1  rosmfnoc 
-rw-rw-r—  1  rosmfnoc 
-rw-rw-r-  1  rosmfnoc 
-rw-rw-r—  1  rosmfnoc 
-rw-rw-r—  1  rosmfnoc 


528482304 

Sep  30 

124583936 

Sep  30 

23724032 

Sep  30 

8760 

Sep  30 

1992  c3ave0885.tar 
1992  c3grid0885.tar 
1992  c3spec0885.tar 


1992  hccn32 
5373952  Sep  30  1992  lsifil0885.tar 

172097536  Sep  30  1992  shist0885.tar 

124464  Sep  30  1992  trpfil 


The  complete  NOGAPS  AMIP  history  data  base  is  located  on  the  SSC  YMP8128 
(designated  LSC)  in  the  master  directory  '/u/b/rosmond/pcmdi'.  In  this  directory  are  subdirectories 
containing  the  120  monthly  histories  from  January  1979  through  December  1988.  The  naming 
convention  of  the  subdirectories  is  'pcmdiMMYY',  so  for  example  'pcmdi0885‘  contains  the 
history  for  August  1985.  Table  2  displays  the  contents  of  'pcmdi0885'  using  the  UNIX  command 
'Is  -la'.  Note  that  'tar'  files  are  identified  with  a  ’.tar'  extension.  The  size  of  each  file  is  given 
in  bytes;  all  files  for  a  month  total  about  800  megabytes.  Table  3  shows,  using  the  UNIX 
command  'Is',  an  example  of  the  expanded  contents  of  the  tar  files.  Except  for  the  two  small 


5 


Table  3.  Expanded  contents  of  tar  files  in  AMIP  directory  "pcmdi0885". 


c3ave057696 

c3grid057966 

c3grid058428 

c3spec058140 

Is  if 8057846 

Is  88058308 

shist058020 

c3ave057720 

c3grid057972 

c3grid058434 

c3spec058146 

Isif8057852 

Is 8805831 4 

shist058026 

c3av»0S7744 

c3grid057978 

c3grid058440 

c3spec058152 

Is  if 8057858 

IS88058320 

shist058032 

c3ave057768 

c3grid057984 

c3spec057696 

c3spec058158 

Isif8057864 

Is 88058326 

shist058038 

c3ave057792 

c3grid057990 

c3spec057702 

c3spec058164 

Isif8057870 

IS88058332 

shist058044 

c3ave057816 

c3grid057996 

c3spec057708 

c3spec058170 

Isif8057876 

Is 88058338 

shist058050 

c3avo057840 

c3grid058002 

c3spec057714 

c3spec0S8176 

Isif 8057882 

Is88058344 

shist058056 

c3ave057864 

c3grid058008 

c3spec057720 

c3spec058182 

Is  if 8057888 

Is88058350 

shist058062 

c3ave057888 

c3grid058014 

c3spec057726 

c3spec058188 

Isif8057894 

Is  if 8058356 

shist058068 

c3ave057912 

c3gnd058020 

c3spec057732 

c3spec058194 

Isif 8057900 

Is88058362 

shist058074 

c3ave057936 

c3grid058026 

c3spec057738 

c3spec058200 

Isif8057906 

Is88058368 

shtst058080 

c3ave057960 

c3grid058032 

c3spec057744 

c3spec058206 

lsif805791 2 

Is88058374 

shist058086 

c3ave057984 

c3grid058038 

c3spec057750 

C3spec058212 

Is  if 8057918 

Is 88058380 

shtst058092 

c3ave058008 

c3grid058044 

c3spec057756 

c3spec058218 

Isif 8057924 

IS88058386 

shist058098 

c3ave058032 

c3grid058050 

c3spoc057762 

c3spec058224 

Is  if 8057930 

Is 88058392 

shist058104 

c3ave058056 

c3grid058056 

c3spec057768 

c3spec058230 

Isifit057936 

Is88058398 

shist058110 

c3ave058080 

c3grid058062 

c3spec057774 

c3spec058236 

Is  if 8057942 

IS88058404 

shist058116 

c3ave058104 

c3gnd058068 

c3spee057780 

c3spec058242 

Is  if 8057948 

Is8805841 0 

shtst058122 

c3ave058128 

c3grid058074 

c3spec057786 

c3spec058248 

Isif8057954 

IS8B058416 

shist058128 

c3ave058152 

c3gnd058080 

c3spec057792 

c3spec058254 

Isif 8057960 

ls88058422 

shist058134 

c3ave058176 

c3g  ncl058086 

c3spec057798 

c3spec058260 

Is  if 8057966 

Is 88058428 

shist058140 

c3ave058200 

c3gnd058092 

c3spec057804 

c3spec058266 

Isif 8057972 

IS88058434 

shist058146 

c3ave058224 

c3grid058098 

c3spec057810 

c3spec058272 

Isif 8057978 

Isif  B058440 

shist058152 

c3ave058248 

c3grid0581 04 

c3spec057816 

c3spec058278 

Isif8057984 

shist057696 

shist058158 

c3ave058272 

c3grid058110 

c3spec057822 

c3spec058284 

Isif 8057990 

shist057702 

shist058164 

c3ave058296 

c3grid058116 

c3spec057828 

c3spec058290 

Is  if 8057996 

shist057708 

shist058l70 

c3ave058320 

c3grid058122 

c3spec057834 

c3spec058296 

lsiT8058002 

shist057714 

shist058176 

c3ave058344 

c3grid058128 

c3spec057840 

c3spec058302 

Isif 8058008 

shist057720 

shist058182 

c3ave058368 

c3grid058134 

c3spec057846 

c3spec058308 

Isif805801 4 

shist057726 

shist058188 

c3ave058392 

c3grid058140 

c3spec057852 

c3spec058314 

Isif 8058020 

shist057732 

shist058194 

c3ave058416 

c3gnd058146 

c3spec057858 

c3spec058320 

Isif 8058026 

shist057738 

shist058200 

c3ave058440 

c3grid058152 

c3spec057864 

c3spec058326 

Isif8058032 

shist057744 

shist058206 

c3grid057696 

c3grid058158 

c3spec057870 

c3spec058332 

Is  if 8058038 

shist057750 

shist058212 

c3gnd057702 

c3grid058164 

c3spec057876 

c3spec058338 

Isif 8058044 

shist057756 

shist058218 

c3grid057708 

c3grid058170 

c3spec057882 

c3spec058344 

Isif 8058050 

shist057762 

shis(058224 

c3grid057714 

c3grid058176 

c3spec057888 

c3spec058350 

Isif8058056 

shist057768 

shist058230 

c3grid057720 

c3grid058182 

c3spec057894 

c3spec058356 

Isif 8058062 

shist057774 

shist058236 

c3grid057726 

c3grid058188 

c3spec057900 

c3spec058362 

Is  if  il058068 

shist057780 

shist058242 

c3grid057732 

c3grid058194 

c3spec057906 

c3spec058368 

Isif8058074 

shist057786 

shist058248 

c3grid057738 

c3grid058200 

c3spec057912 

c3spec058374 

Isif8058080 

shist057792 

shist058254 

C3grid057744 

c3grid058206 

c3spec05791 8 

c3spec058380 

Isif8058086 

shist057798 

shist058260 

c3grid057750 

c3grid058212 

c3spec057924 

c3spec058386 

Isif 8058092 

shisf057804 

shist058266 

c3gnd057756 

c3grid058218 

c3spec057930 

c3spec058392 

Isif 8058098 

shist057810 

shist058272 

c3grid057762 

c3grid058224 

c3spec057936 

c3spec058398 

Isif 80581 04 

shist057816 

shist058278 

c3grid057768 

c3grid058230 

c3spec057942 

c3spec058404 

Isif 80581 10 

shist057822 

shist058284 

C3grid057774 

c3grid058236 

c3spec057948 

c3spec058410 

Isifil058116 

shisf057828 

shist058290 

c3grid057780 

c3grid058242 

c3spec057954 

c3spec058416 

Isif8058122 

shist057834 

shist058296 

c3grid057786 

c3grid058248 

c3spec057960 

c3spec058422 

lsif8058 128 

shist057840 

shist058302 

c3grid057792 

c3grid058254 

c3spec057966 

c3spec058428 

Isif8058134 

shist057846 

shist058308 

c3grid057798 

c3grid058260 

c3spec057972 

c3spec058434 

Isif 80581 40 

shist057852 

shist058314 

c3grid057804 

c3grid058266 

c3spec057978 

c3spec058440 

Is  if 80581 46 

shist057858 

shist058320 

c3grid057810 

c3grid058272 

c3spec057984 

hccn32 

Isif8058152 

shist057864 

shist058326 

c3grid057816 

c3grid058278 

c3spec057990 

Isifit057696 

Isifil058158 

shist057870 

shist058332 

c3grid057822 

c3grid058284 

c3spec057996 

Isif 1057702 

Isif8058164 

shist057876 

shist058338 

c3grid057828 

c3grid058290 

c3spec058002 

lsif«057708 

Isif8058170 

shist057882 

shist058344 

c3grid057834 

c3grid058296 

c3spec058008 

Isifil05771 4 

Isif 8058 176 

shist057888 

shist058350 

c3grid057840 

c3grid058302 

c3spec05801 4 

Isifil057720 

Isrf80581 82 

shist057894 

shist058356 

c3gnd05784€ 

c3grid058308 

c3spec058020 

Isif 1057726 

Is  if 8058188 

shist057900 

shist058362 

c3grid057852 

c3grid058314 

c3spec058026 

Isifil057732 

lsifB058194 

shist057906 

shist058368 

c3grid057858 

c3grid058320 

c3spec058032 

Is  if  W5 7738 

Isif8058200 

shist057912 

shist058374 

c3grid057864 

c3grid058326 

c3spec058038 

Isif 10577 44 

Isif 8058206 

shist057918 

shist058380 

c3grid057870 

c3grid058332 

c3spec058044 

Isif  KJ57750 

Isifil058212 

shist057924 

shist058386 

c3grid057876 

c3grid058338 

c3spec058050 

Isifil057756 

Isif8058218 

shist057930 

shist058392 

c3grid057882 

c3grid058344 

c3spec058056 

Is  if  8057762 

Isif 8058224 

shist057936 

shist058398 

c3grid057888 

c3g  rid058350 

c3spec058062 

Isify057768 

Isif 8058230 

shist057942 

shist058404 

c3grid057894 

c3grid058356 

c3spec058068 

Is  if  8057774 

Is  if 8058236 

shist057948 

shist058410 

c3grid057900 

c3grid058362 

c3spec058074 

Isif 8057780 

Isif 8058242 

shist057954 

shist058416 

c3grid057906 

c3grid058368 

c3spec058080 

Isif 8057786 

Isif 8058248 

shist057960 

shist058422 

c3grid057912 

c3grid058374 

c3spec058086 

Isif 8057792 

Isif8058254 

shist057966 

shist058428 

c3grid057918 

c3grid058380 

c3spec058092 

lsKU057798 

Isif 8058260 

shist057972 

shist058434 

c3grid057924 

c3grid058386 

c3spec058098 

ls#8057804 

Isif8058  2  66 

shist057978 

shist058440 

c3grid057930 

c3grid058392 

c3spec058104 

Is  if 805781  0 

lsif8058272 

shist057984 

trpf8 

c3grid057936 

c3grid058398 

c3spec058110 

Isif8057816 

Is  if 8058278 

shist057990 

C3grid057942 

c3grid058404 

c3 spec0581 16 

Isif8057822 

lsiT8058284 

shist057996 

c3grid057948 

c3gnd058410 

c3spec058122 

Is  if 8057828 

Isif8058290 

shist058002 

c3gnd057954 

c3grid058416 

c3spec058128 

Isif8057834 

Isif8058296 

shtsl058008 

c3grid057960 

c3grid058422 

c3spec058134 

Isif 8057840 

Isif 8058302 

shist05801 4 
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auxiliary  files  'hccn32'  and  'trpfil',  all  the  files  have  a  integer  counter  appended  to  the  end  of  the 
text  part  of  the  filenames.  This  number  is  the  number  of  hours  (tau)  since  the  beginning  of  the 
AMIP  integration  on  January  1,  1979  00  UTC.  The  software  described  in  the  appendix  contains 
Fortran  subroutines  for  converting  between  date-time-group  and  tau.  For  the  6  hour  interval 
output  files,  the  data  contained  is  a  snapshot  of  the  model  state  at  this  forecast  time.  The  24  hour 
interval  files  (c3ave...)  contain  daily  means  of  various  model  parameters  averaged  over  every 
model  time  step  for  the  24  hour  period  ending  at  the  tau  value  in  the  filename.  Each  file's 
contents  are  described  in  detail  in  the  next  section. 

4.  NOGAPS  HISTORY  FILE  DESCRIPTIONS 


The  NOGAPS  AMIP  model  history  files  contain  a  wide  variety  of  output  variables.  When 
NOGAPS  runs  as  a  NWP  model  the  model  history  output  is  in  the  form  of  instantaneous 
snapshots  of  the  model  atmosphere  state  every  6  hours.  The  NWP  output  is  contained  in  the  files 

with  names  beginning  with  the  text:  'shist....',  'c3spec....\  'c3grid....',  and  'lsifil _ with  appropriate 

'tau'  values  appended  as  shown  in  table  3.  For  the  AMIP  experiment  the  NWP  history  was 
augmented  with  daily  means  and  standard  deviations  of  several  model  variables  and  parameters. 
These  data  are  contained  in  the  'c3ave....'  files  shown  in  Table  3.  Because  the  fundamental  I/O 
mechanism  used  for  these  files  is  Fortran  direct-access,  all  records  within  a  file  must  be  the  same 
length.  The  files  'shist....'  and  'c3spec....'  contain  NOGAPS  history  data  that  are  stored  as 
spherical  harmonic  coefficients  and  have  record  lengths  of  25936  bytes.  The  files  ’c3grid....', 
'lsifil....',  and  'c3ave....'  contain  NOGAPS  history  data  on  the  NOGAPS  Gaussian  transform  grid 
and  have  record  lengths  of  115216  bytes.  Tables  4,  5,  6,  7,  and  8  give  the  itemized  contents  of 
'shist....',  'c3spec....',  'c3grid....',  'lsifil....',  and  'c3ave....'. 


Table  4.  File  'shist....'  contains  spherical  harmonic  coefficient  data  for  two  consecutive 
time  levels  at  the  model  history  output  time.  Record  lengths  are  25936  bytes. 


'shist'  NOGAPS  Variable 

record  no(s). 

Vorticity  at  current  time  step 

1-18 

Divergence  at  current  time  step 

19-36 

Temperature  at  current  time  step 

37-54 

Specific  humidity  at  current  time  step 

55-72 

Terrain  pressure  at  current  time  step 

73 

Vorticity  at  previous  time  step 

74-91 

Divergence  at  previous  time  step 

92-109 

Temperature  at  previous  time  step 

110-127 

Specific  humidity  at  previous  time  step 

128-145 

Terrain  pressure  at  previous  time  step 

146 
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Table  5.  File  'c3spec....'  contains  auxiliary  spherical  harmonic  coefficient  data  used 
internally  by  NOGAPS.  Record  lengths  are  25936  bytes. 


'c3spec'  NOGAPS  Variable 

record  no(s). 

Radiation  heating  rate 

1-18 

Laplacian  of  terrain  geopotential 

19 

Terrain  pressure  tendency 

20 

Table  6.  File  'c3grid....'  contains  2-dimensional  Gaussian  grid  fields  of  parameters  describing 
characteristics  of  NOGAPS  surface  features  and  physical  processes.  Record  length 
is  115216  bytes. 


'c3grid'  NOGAPS  Variable 

record  no(s). 

Terrain  geopotential  (m/sec)*  *2 

1 

Standard  deviation  of  terrain  geopotential  (m/sec)* *2 

2 

Ground  temperature  (K) 

3 

Ground  wetness  climotology  (0-1) 

4 

Ground  temperature  climotology  (K) 

5 

Ground  wetness  (0-1) 

6 

Surface  albedo  (0-1) 

7 

Surface  rougness  (m) 

8 

Snow  depth  (mm  h2o) 

9 

Surface  solar  heat  flux  (W/m**2) 

10 

Surface  longwave  heat  flux  (W/m**2) 

11 

Surface  evaporation  (kg/m* *2) 

12 

Surface  moist  static  energy  flux  (W/m**2) 

13 

Cumulus  cloud  base  mass  flux  (kg/m** 2/sec) 

14 

Surface  drag  (Nts/m**2) 

15 

Accumulated  convective  precip  (mm) 

16 

Accumulated  stable  precip  (mm) 

17 

Total  precipitation  rate  (mm/sec) 

18 

Top  of  atmosphere  net  solar  flux  (W/m**2) 

19 

Top  of  atmosphere  longwave  flux  (W/m**2) 

20 

Lifting  condensation  level  (mb  above  surf) 

21 

Top  of  deepest  convective  cloud  (mb) 

22 

Cumulus  cloud  fraction  (0-1) 

23 

Unused 

24 

Table  7.  File  Tsifil _ '  contains  the  NOGAPS  land/sea/ice  table. 

Record  length  is  115216  bytes. 


'lsifil'  NOGAPS  Variable 

record  no(s). 

Opensea/barel  and/seaice/landice  (0 : 1 :2 : 3) 

1 

Besides  the  tau-tagged  history  files,  Table  3  shows  two  small  auxiliary  files,  'hccn32'  and 
'trpfir,  which  do  not  have  a  'tau'  value  appended.  File  'hccn32'  contain  the  set  of  model 
configuration  parameters  (e.g.,  distribution  of  vertical  layers)  and  physical  constants  (e.g., 
acceleration  of  gravity)  which  were  used  for  the  AMIP  integrations.  During  all  model  restarts 
this  file  is  read  to  ensure  that  these  quantities  are  unchanged  for  the  entire  integration  period. 

File  'trpfil'  contains  three  Gaussian  grid  fields,  defined  at  the  beginning  of  the  forecast, 
used  to  do  the  reduction  to  sea  level  necessary  to  output  fields  on  constant  pressure  surfaces. 
The  first  field  is  the  'pdiff  array,  which  is  the  pressure  difference  between  sea  level  and  terrain 
level  at  the  initial  forecast  time,  which  for  the  AMIP  integration  was  1  Jan  1979  00  UTC. 
'Pdiff  ranges  from  zero  millibars  over  the  oceans  to  about  500  millibars  over  the  Himalayas.  The 
second  field  is  a  surface  to  100  millibar  layer  mean  atmospheric  temperature.  Changes  in  this 
deep  layer  mean  during  the  forecast  are  used  to  adjust  the  'pdiff  values  to  allow  for  the  effect 
of  warming  or  cooling,  on  the  fictitious  subterranean  air  column  between  sea  level  and  terrain 
level.  The  third  field  is  the  initial  1000  mb  temperature.  This  field,  combined  with  the  second 
field  above,  are  used  in  short  term  NWP  forecasts  to  accurately  predict  sea  level  pressure 
changes. 

Users  of  the  NOGAPS  model  history  will  probably  have  their  own  ideas  about  how  to  do 
the  reduction  to  sea  level  process.  However,  for  the  purposes  of  looking  at  long  term  means, 
we  suggest  that  the  'pdiff  array  above  be  added  to  time  averaged  terrain  pressures  without 
worrying  about  the  subtleties  of  temperature  correcting  used  in  short  term  NWP  forecasts.  Sea 
level  pressure  prediction  under  high  terrain  is  certainly  not  a  priority  question  in  long  term 
atmospheric  simulation  and  should  be  kept  as  simple  as  possible. 

5.  SUMMARY 

Anyone  with  an  account  on  the  SSC  computer  system  can  access  the  NOGAPS  AMIP 
history  data  base.  All  files  have  read  access  for  everyone.  There  is  no  anonymous  FTP  access 
to  SSC,  so  unfortunately  outside  access  is  not  available.  If  non-SSC  parties  are  interested  in  the 
AMIP  data,  other  arrangements  may  be  possible,  however. 
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Table  8.  File  'c3ave....'  contains  Gaussian  grid  fields  of  daily  means  and  variances  of 
several  NOGAPS  predicted  quantities.  Record  length  is  115216  bytes. 


'c3ave'  NOGAPS  Variable 

record  no(s). 

Latent  heat  flux  (W/m**2) 

l 

Sensible  heat  flux  (W/m**2) 

2 

Top  of  atmosphere  solar  flux  (W/m**2) 

3 

Top  of  clear  atmosphere  solar  flux  (W/m**2) 

4 

Surface  solar  flux  (W/m**2) 

5 

Surface  solar  flux,  clear  atm  (W/m**2) 

6 

U-comp,  surface  drag  (Nts/m**2) 

7 

V-comp,  surface  drag  (nts/m**2) 

8 

Top  of  atmosphere  longwave  flux  (W/m**2) 

9 

Surface  longwave  flux  (W/m**2) 

10 

Top  of  clear  atmosphere  longwave  flux  (W/m**2) 

11 

Surface  longwave  flux,  clear  atm  (W/m**2) 

12 

Surface  air  temperature  (K) 

13 

Spectral  truncation  surface  moisture  flux  (W/m**2) 

14 

Convective  precipitation  (mm) 

15 

Stable  precipitation  (mm) 

16 

Terrain  pressure  (mb) 

17 

Variance  of  surface  air  temperature  (K**2) 

18 

Ground  temperature  (K) 

19 

Variance  of  ground  temperature  (K**2) 

20 

Diabatic  heating  (K/day) 

21-38 

Diabatic  moistening  (kg/kg/day) 

39-56 

Diabatic  U-comp  change  (m/sec/day) 

57-74 

Diabatic  V-comp  change  (m/sec/day) 

75-92 

Precipitation  heating  (K/day) 

93-110 

Precipitation  moistening  (kg/kg/day) 

111-128 

Type  2  cloud  forcing  stencil 

129 

Cumulus  friction,  U-comp  (m/sec/day) 

130-146 

Cumulus  friction,  V-comp  (m/sec/day) 

147-164 

Total  cloud  fraction  (0-1) 

165 

Total  precipitable  water  (kg/m**2) 

166 

Cumulus  ensemble  cloud  base  mass  flux  (kg/m**2/day) 

167-182 

Gravity  wave  drag,  U-comp  (Nts/m**2) 

183-200 

Gravity  wave  drag,  V-comp  (Nts/m**2) 

201-218 

Solar  heating  (K/day) 

219-236 

Longwave  heating  (K/day) 

237-254 

Solar  heating,  clear  atmos  (K/day) 

255-272 

Longwave  heating,  clear  atmos  (K/day) 

273-290 

Stable  cloud  fraction  (0-1) 

291-308 

Convective  cloud  fraction  (0-1) 

309-326 

Temperature  change,  horizontal  diffusion  (K/day) 

327-344 

Moisture  change,  horizontal  diffusion  (kg/kg/day) 

345-362 

U-comp  change,  horizontal  diffusion  (m/sec/day) 

363-380 

V-comp  change,  horizontal  diffusion  (m/sec/day) 

381-398 
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In  the  following  appendix  several  subroutines  for  I/O  and  data  management  and 
manipulation  are  listed.  Also  listed  is  a  program  used  to  process  the  data  from  one  of  the  history 
directories  into  a  standard  set  of  monthly  mean  quantities  requested  by  PCMDI  from  each  of 
the  AMIP  participants.  Users  should  find  it  a  useful  example  to  modify  or  base  their  own 
AMIP  data  processing  program.  All  source  code  of  the  listed  software  is  in  the  directory 
'a/rosmond/pcmdi/sources'.  Users  are  warned  that  this  software  is  written  for  CRAY  systems  and 
has  some  unique  requirements  unlikely  to  be  available  on  other  platforms.  Specifically,  there  are 
several  references  to  CRAY  scientific  library  routines  such  as  FFTs,  and  the  non-standard 
dynamic  array  extension  is  used  in  several  subroutines.  If  necessary,  however,  the  author  is 
willing  to  provide  some  advice  about  software  conversion. 
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APPENDIX  -  SOFTWARE  LISTINGS 


subroutine  cons( rad, radsq, msort, l sort, ml  sort ,eps4,ci m, weight, si nl 
*,  ocos,poly,dpoly) 
c 

c  subroutine  to  create  a  nunber  of  essential  constant  operator  arrays 
c  for  NOGAPS  history  file  processing  (T47) 
c 

c  output : 
c 

c  rad:  radius  of  earth:  6371000.0 
c  radsq:  rad* rad 

c  msort:  zonal  wavenumber (m)=  msortfspherical  harmonic  index(ml)) 
c  Isort:  total  wavenumber( l )=  lsort(spherical  harmonic  index(ml)) 
c  mlsort:  spherical  harmonic  index(ml)=  mlsort(m,l) 
c  eps4:  spherical  harmonic  laplacian  operator 
c  cim:  fourier  coefficient  zonal  wavenumber  operator 
c  weight:  gaussian  quadrature  weights 
c  sinl:  sine  of  gaussian  latitudes 
c  ocos:  1.0/cos(latitude)**2 
c  poly:  associated  legendre  polynomials 
c  dpoly:  d(poly)/d(sinl ) 
c 

c  ****************************************************************** 
c 

parameter  (im=144,  lm=18,  lpout=6) 
parameter  (jm=  im/2,  im2=  im/2,  jm2=  jm/2) 

parameter  (]trun=  2*((1+( im- 1 )/3)/2),  mlmax=  ( jtrun+1 )*jtrun/2) 
parameter  (imlm=im*lm,  imjm=im*jm,  idim2=  mlmax*2) 
c 

common/  fftcom/  trigs(im,2), ifax(19) 
c 

dimension  msort (mlmax), lsort(mlmax),m(sort( jt run, jtrun) 

*,  eps4(mlmax),cim(mlmax),weight( jm),sinl( jm),ocos{ jm) 

*,  polylmlmax, jm2),dpoly(mlmax, jm2) 
c 

c  generate  sorting  arrays 
c 

cal l  sortml( jtrun, mlmax, msort, Isort, ml  sort) 
c 

c  generate  laplacian  operator  and  fourier  coefficient  zonal 
c  wavenumber  coefficient 
c 

rad=  6371000.0 
radsq=  rad* rad 
do  150  ml=1, mlmax 
rl=  Isort(ml) 
rm=  msort(ml)-1 
rlm=  rl-1.0 

if(msort(ml).eq.1)  rm=  0.0 
if(lsort(ml).eq.1)  rlm=  0.0 
eps4(ml)=  rl*rlm/radsq 
cim(ml)=  rm 
150  continue 
c 

c  generate  gaussian  quadrature  weights  and  latitudes 
c 

call  gausl3( jm, -1 .0, 1 .0, weigh t, sinl ) 
c 

do  155  j=1,jm 

ocos(j)=  1 .0/(1 ,0-sinK  j)*sinl(  j ) ) 

155  continue 
c 

c  generate  legendre  functions 


A-  1 


c 

call  lgndr<  jm2,  jtrun.mlmax.uilsort.poly.dpoly.sinl) 
c 

c  generate  CRAY  library  fft  coefficients 
c 

call  fftfaxd'm, ifax, trigs) 
c 


return 

end 


A-2 


rilHt 


subroutine  daynum(cidtg,my,rmi,md,mh,mdy,mhy) 
c 

c  subroutine  to  convert  a  date-time-group  to  equivalent 
c  year/month/day/hour  set 
c 

c  input: 

c 

c  cidtg:  ascii  date-time-group  (YYMMDOHH) 
c 

c  output : 

last  two  digits  of  year 
month  of  year 
day  of  month 
hour  of  day 

:  julian  date  (sideral  day) 

:  hour  of  month 


£*********************************************************** 

c 

character*8  cidtg 
c 

dimension  month(12,2) 

data  month/0, 31, 59, 90, 120, 151, 181 ,212, 243,273, 304, 334, 
*  0,31,60,91,121,152,182,213,244,274,305,335/ 

c 

read(cidtg,1)  my,mm,md,mh 
1  format(4i2) 

j  =  1 

if(mod(my,4).eq.O)  j  =2 
mdy  =  month(mm, j )+md 
mhy  =  (mdy-1 )*24+mh 
return 
end 
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subroutine  dtgf ix(cidtg, cndtg, ndif ) 
c 

c  subroutine  to  compute  new  date-time-group  (dtg)  from  base 
c  dtg  and  offset  period  in  hours 
c 

c  i nput : 

c 

c  cidtg:  input  dtg  (YYMMODHH) 
c  ndif:  offset  time  in  hours 
c 

c  output : 

c 

c  cndtg:  output  dtg 
c 

c**«*************************************************»***** 
characters  cidtg,  cndtg 
dimension  myc(4) 
data  myc/8784, 3*8760/ 
call  daynun(cidtg,my,nin,md,mh,mdy,mhy) 
c 

c  increment  (decrement)  hour  of  the  year  (mhy)  by  ndif, 
c  test  for  change  of  year,  adjust  if  necessary,  reform  to  ndtg 
c 

if(my.eq.O)  my=  100 
newmhy  =  mhy+ndif 
if(newmhy.ge.O)  go  to  10 
1  my  =  my- 1 

idx  =  mod(my,4)+1 
newmhy  =  myc( idx)+newmhy 
if (newmhy. ge.O)  go  to  20 
go  to  1 
c 

10  idd  =  mod(my,4)*1 

if (newmhy. It. myc(idd))  go  to  20 
newmhy  =  newmhy-myc( idd) 
my=  my+1 
go  to  10 
c 

20  my=  mod(my,100) 

call  noday( cndtg, my, newmhy) 

return 

end 
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subroutine  gausl3  (n,xa,xb,wt,ab) 
c 

c  subroutine  to  compute  gaussian  latitudes  and  quadrature  weights 
c  for  gaussian  quadrature  integrations 
c 

c  **************************************************************** 
c 

c  weights  and  abscissas  tor  nth  order  gaussian  quadrature  on  (xa,xb). 
c  input  arguments 

c  n  -the  order  desired  (number  of  gaussian  latitudes) 

c  xa  -the  left  endpoint  of  the  interval  of  integration  (-1.0  for  spole) 

c  xb  -the  right  endpoint  of  the  interval  of  integration(1 .0  for  npole) 

c  output  arguments 

c  ab  -the  n  calculated  abscissas 

c  wt  -the  n  calculated  weights 

c 

c  ****************************************************************** 
c 

implicit  double  precision  (a-h,o-z) 
c 

real  ab(n)  ,wt(n) 

c 

c  machine  dependent  constants — 

c  tol  -  convergence  criterion  for  double  precision  iteration 


c 

pi 

-  given  to  15  significant 

digi ts 

c 

cl 

-  1/8 

these  are  coefficients  in  mcmahon"s 

c 

c2 

-  -31/(2*3*8**2) 

expansions  of  the  kth  zero  of  the 

c 

c3 

-  3779/(2*3*5*8**3) 

bessel  function  j0(x)  (cf.  abramowitz 

c 

c4 

-  -6277237/(3*5*7*8**5) 

handbook  of  mathematical  functions). 

c 

u 

-  ( 1 - ( 2/pi )**2)/4 

c 

data  tol/1 .d-15/,pi/3.141 59265358979/ , u/ . 1 48678816357662/ 
data  c 1 , c2 , c3 , c4/ . 1 25 , - . 080729 1 66666667, . 246028645833333 , 

1  -1.82443876720609  / 

c 

c  maximum  nunber  of  iterations  before  giving  up  on  convergence 
c 

data  maxit  /5/ 
c 

c  arithmetic  statement  function  for  converting  integer  to  double 
c 

dbli(i)  =  dble(float(i)) 
c 

ddif  =  .5d0*(dble(xb)-dble(xa)) 
dsum  =  .5d0*(dble(xb)+dble(xa)) 
if  (n  .gt.  1)  go  to  101 
ab(1)  =  0. 
wt(1)  =  2.*ddif 
go  to  107 
101  continue 

nnpl  =  n*(n+1 ) 

cond  =  1./sqrt((.5+float(n))**2+u) 
l im  =  n/2 
c 

do  105  k=1, lim 

b  =  (float(k)-.25)*pi 
bisq  =  1./(b*b) 
c 

c  rootbf  approximates  the  kth  zero  of  the  bessel  function  j0(x) 
c 

rootbf  =  b*(1 .+bisq*(c1*bisq*(c2+bisq*(c3+bisq*c4)))) 
c 

c  initial  guess  for  kth  root  of  tegendre  poly  p-sub-n(x) 

c 

dzero  =  cos( rootbf *cond) 
do  103  i=1, maxit 
c 
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dpm2  =  l.dO 
dpml  *  dzero 


recursion  relation  for  legendre  polynomials 
do  102  nn=2,n 

dp  =  <dbli(2*nn-1)*dzero*dpm1-dbli(nn-1)*dpm2)/dbti(nn) 
dpm2  =  dpml 
dpml  =  dp 

102  continue 

dtmp  *  1 .d0/(1 .d0-dzero*dzero) 

dppr  =  dbli(n)*(dpm2-dzero*dp)*dtmp 

dp2pri  *  (2.d0*dzero*dppr-dbl i (nnpl )*dp)*dtmp 

drat  =  dp/dppr 

cubical ly-convergent  iterative  improvement  of  root 

dzeri  *  dzero-drat*(1.d0+drat*dp2pri/(2.d0*dppr)) 
ddun=  dabs  (dzeri -dzero) 
if  (ddum  . le.  tol)  go  to  104 
dzero  *  dzeri 

103  continue 
print  504 

504  formatOx, 1  in  gausl3,  convergence  failed1) 

104  continue 

ddifx  =  ddif*dzero 
ab(k)  =  dsun-ddifx 

nt(k)  *  2.dO*(1.dO-dzero*dzero)/(dbli(n)*dpm2)**2*ddif 
i  =  n-k+1 

ab(i)  =  dsum-ddifx 
wt(i)  =  wt(k) 

105  continue 

if  (mod(n,2)  .eq.  0)  go  to  107 

ab(lim»1)  =  dsum 

nml  *  n-1 

dprod  =  n 

do  106  k=1,nm1,2 

dprod  =  dbl i (nml -k)*dprod/dbl i (n-k) 

106  continue 

wt(lim+1)  =  2.d0/dprod**2*ddi f 

107  return 
end 
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subroutine  lgrtdr( jm2, jtrun, mlmax, ml sort, poly ,dpoty,s ini) 
c 

c  **  neprf  ***,  programmer  tom  rosmond,  august  3,  1987 
c 

c  generate  legendre  polynomials  and  their  derivatives  on  the 
c  gauss i an  latitudes 
c 

c  ref=  belousov,  s.  1.,  1962=  tables  of  normalized  associated 
c  legendre  polynomials,  pergamon  press,  new  york 

c 

c  ********************************************************** 

c 

c  input: 

c 

c  jm2:  number  of  gaussian  latitudes  between  spole  and  equator 
c  jtrun:  zonal  wavenumbers  (48  for  T47  model) 
c  mlmax:  number  of  degrees  of  freedom  in  spherical  harmonic 
c  variables  =  jtrun*( jtrun+1 ) 

c  mlsort:  sorting  array  for  locating  spherical  harmonics  as  a 
c  function  of  zonal  waveno(m)  and  total  waveno(l) 

c  sinl:  sine  of  gaussian  latitudes 
c 

c  output : 

c 

c  poly:  associated  legendre  polynomials 
c  dpoly:  d(poly)/d(sinl ) 
c 

c  ************************************************************ 
c 

dimension  poly(mlmax, jm2) , dpoly (mlmax, jm2),sinl( jm2) 

*,  mlsort( jtrun, jtrun) 
c 

dimension  pnm(49,48) ,dpnm(49,48) 
c 

c  sinl  is  sinClatitude)  =  cos(colatitude) 

c  prm(np.mp)  is  legendre  polynomial  p(n,m)  with  np=n+1,  mp=m+1 
c  pnm(mp,np+1 )  is  x  derivative  of  p(n,m)  with  np=n+1,  mp=m+1 
c 

jtrunp=  jtrun+1 
do  1001  j=1 , jm2 
xx=  sinl(j) 
sn=  sqrt(1 -0-xx*xx) 
sn2i  =  1 .0/(1 .0  -  xx*xx) 
rt2=  sqrt(2.0) 
cl  =  rt2 
c 

pnm(1,1)  =  1.0/rt2 

theta=-atan(xx/sqrt(1 ,0-xx*xx))+2.0*atan( 1.0) 
c 

do  20  n=1, jtrun 
np  =  n  ♦  1 
fn=n 

fn2  =  fn  +  fn 
fn2s  =  fn2*fn2 

c  eq  22 

c1=  c1*sqrt(1.0-1.0/fn2s) 
c3=  c1/sqrt(fn*(fn+1.0)) 
ang  =  fn*theta 
si  =  0.0 
s2  =  0.0 
c4  =  1.0 
c5  =  fn 
a  =  -1.0 
b  =  0.0 
c 

do  27  kp=1,np,2 
k  *  kp  -  1 
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s2=  s2+c5*sin(ang)*c4 
if  (k.eq.n)  c4  =  0.5*c4 
s1=  s1+c4*cos(ang) 
a  =  a  ♦  2.0 
b  =  b  ♦  1.0 
fk=k 

ang  =  theta*(fn  -  fk  -  2.0) 
c4  *  (a*(fn  -  b  +  1.0)/(b*(fn2  -  a)))*c4 
c5  =  c5  -  2.0 
27  continue 
c  eq  19 

pnmCnp, 1 )  =  s1*c1 

c  eq  21 

pnm(np,2)  =  s2*c3 
20  continue 
c 

do  4  mp=3, jtrunp 
m  =  mp  -  1 
fm=  m 

fml  =  fm  -  1.0 
fm2  =  fm  -  2.0 
fm3  =  fm  -  3.0 
c6=  sqrt(1 .0+1 .0/(fm+fm)) 
c  eq  23 

pnm(mp,mp)  a  c6*sn*pnm(m,m) 
if  (mp  -  jtrunp)  3,4,4 

3  continue 

nps  =  mp  ♦  1 
c 

do  41  np=nps, jtrunp 
n  =  np  -  1 
fn=  n 

fn2  =  fn  +  fn 

c7  =  (fn2  ♦  1.0)/(fn2  -  1.0) 
cfl  =  (fml  ♦  fn)/((fm  +  fn)*(fm2  ♦  fn)) 
c=  sqrt((fn2+1.0)*c8*(fm3+fn)/(fn2-3.0)) 
d=  - sqr t ( c7*c8* ( f n- fml ) ) 
e=  sqrt(c7*(fn-fm)/(fn+fm)) 
c  eq  17 

prm(np,mp)  a  c*pnm(np-2,mp-2) 

1  ♦  xx*(d*pnm(np-1 ,mp-2)  ♦  e*pnm(np  -  1,mp)) 

41  continue 

4  continue 
c 

do  50  mp=1,jtrun 
fm=  mp-1.0 
fms  =  fm*fm 
do  50  np=mp, jtrun 
fnp=  np 

fnp2  =  fnp  +  fnp 

cf  a  (fnp*fnp  -  fms)*(fnp2  -  1.0)/(fnp2  +  1.0) 
cf a  sqrt(cf) 

e  der 

dpnm(np,mp)  =  -sn2i*(cf*pnm(np+1 ,mp)  -  fnp*xx*pnm(np,mp)) 
50  continue 
c 

do  71  m=1, jtrun 
do  71  l=m, jtrun 
ml=  mlsort(m, l ) 
poly(ml,  j)=  pnmd.m) 
dpoly(mt, j)=dpnm(l,m) 

71  continue 

dpoly(1,j)=  0.0 
1001  continue 
return 
end 
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subroutine  nfopen(f i lnam,msg, i read, i tau, itype, ten, iun.fi le, istat) 
c 

c  subroutine  to  open  NOGAPS  history  files 
c 

c  input: 

c 

c  filnam:  path/filename 
c  msg:  flag  to  control  diagnostic  prints 

c  iread:  flag  set  if  opening  with  first  access  to  read  (returns 
c  bad  status  if  no  prexi sting  file) 

c  itau:  forecast  time  in  hours  (if  itau  It  0,  tau  is  irelevant) 
c  itype:  flag  for  32  bit  vs  64  bit  I/O  (itype  .ne.  0:  32  bit) 
c  len:  number  of  data  elements  (words)  to  records 
c 

c  output : 

c 

c  filnam:  input  filnam  with  itau  (6  digits)  appended 
c  istat:  status  (ne  0:  problems  with  open) 
c 

c  ************************************************************** 

c 

character*54  file 
character*48  filnam 
character*8  status 
character*6  ctau 
c 

logical  lex, opn, iread 
c 

lenr=  8*(2+len) 

if(itype.ne.O)  lenr=  8*(2+( len+1 )/2) 
c 

call  chlen(f i Inam, lenc) 
wri te(ctau, 1 (i6.6) ■ )  itau 
c 

if(itau.lt.O)  then 
file®  filnamd: lenc) 
else 

file®  f  i  InamO  :  lenc)//ctau 
endi  f 
c 

inqui re(f i le=f i le,exist=lex,opened=opn) 
if (opn)  return 
c 

if(lex)  then 
status='old‘ 
istat®  0 
else 

status='new' 

c 

if (iread)  then 
istat=2 
return 
endif 
c 

endif 


open(unit=iun,f i le=f i le, access='di rect ' , form®' unformat ted1 
*,  recl=lenr,status®status) 

if(msg.lt.2)  return 

print  100,  status(1:3), filnam, iun.lenr 
100  format(1x,a3, 1x,a54, '  opened  as  unit=',i3,'  :  tenr=',i9,'  bytes  ’) 
return 
end 

subroutine  nfread(fnam, msg, itype, istrt ,nrec, len,x, itau.cdtg, istat) 
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c  subroutine  to  read  NOGAPS  history  files  (but  can  also  be  used  for 
c  general  purpose  10).  file  structure  is  FORTRAN  direct  access,  as  an 
c  option  can  retrieve  data  written  in  IEEE  format  (e.g.  SUN  format)  and 
c  convert  it  to  CRAY  64  bit  data, 
c 

c  formal  parameters 
c 

c  INPUT: 
c 

c  fnam:  character*48  path/filename  to  data  to  be  read,  path  can  be 
c  relative  or  absolute  but  must  be  pre-existing,  actual  filename  can 
c  have  a  indexing  string  appended  depending  of  value  of  formal  parameter 
c  'itau'.  see  example  below, 
c 

c  msg:  flag  to  control  information  prints  about  file  activity 
c  =0,  no  prints 

c  =1,  print  info  about  read  or  write  action 

c  =2,  print  info  about  file  opening,  reading  and  writing  action 

c  =  anything  else,  same  as  2 

c 

c  itype:  data  type 

c  =  0,  do  not  unpack  data,  used  if  file  data  is  already  CRAY  format 

c  =1  (or  -1),  unpack  from  IEEE  integer  to  CRAY  integer 

c  =2  (or  -2),  unpack  from  IEEE  float  pt.  to  CRAY  float  pt. 

c 

c  istrt:  beginning  direct  access  index  of  records  to  be  read  from  fnam 
c 

c  nrec:  number  of  records  (beginning  at  istrt)  to  be  read  from  fnam 
c 

c  len:  number  of  data  items  read  into  array  x.  (note:  since 
c  file  format  is  direct  access,  len  must  be  the  same  for  all  records 
c  on  a  given  fnam) 
c 

c  itau:  counter  index,  if  .GE.  zero  it  is  appended  as  a  character*6 
c  string  to  fnam,  yielding  a  new  filename  which  is  the  actual 
c  one  opened  for  the  requested  data,  see  example  below, 
c 

c  OUTPUT: 
c 

c  x:  array  that  data  is  being  read  into,  should  have  dimension  at  least 
c  (nrec*len) 

c 

c  itau:  counter  index,  integer  item  included  with  data  read  and  returned  to 
c  calling  program,  the  returned  value  of  itau  can  be  compared  to  the  input 
c  value  to  check  data  integrity, 

c 

c  cdtg:  data-time-group  character*8  item  of  data  read  and  returned  to  calling 
c  program, 

c 

c  istat:  returned  status  code, 
c  =0,  no  problems 

c  =1,  bad  read,  data  not  returned 

c  =2  file  missing,  no  data  returned 

c 

c  EXAMPLE: 
c 

c  ******************************************************************** 
c 

!  program  test 

c 

c  program  to  demonstrate  use  of  subroutines  nfread  and  nfwrits 
c  remove  I's  from  col  inn  1  to  get  compilible  code 
c 

!  parameter  (len=  10001) 

!  dimension  x(len,5) 

!  dimension  xfull(len,5) 

!  dimension  ix(len,5) 
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I 

f 

c 

c 

c 

I 

! 

! 

c 

c 

c 


c 

c 

c 

c 

I 

c 

! 

c 

c 

c 

c 

c 

c 

I 

I 

! 

; 

i 

! 

c 

! 

c 

c 

c 

c 

i 

i 

! 

c 

i 

c 

c 

c 

c 

c 

c 

c 

I 

I 

I 

I 


c 

I 

c 


character*48  ifile 
character*8  cdtg 

change  the  following  path/filename  to  suit  your  interests 

if i l e='/a/yourname/datapa th/dummy' 

cdtg*1 92030412' 

n=5 

itau*  10 

generate  some  test  data 

do  1  j=1,n 
do  1  i=1,len 
x(i,j)=  (i+j) 
x(i, j)>-x(if j)**3.888 
xfull(i, j)=  x(i,j> 
ix(t,j)*  <i+j)**2 
1  continue 

this  is  the  largest  integer  preserved  when  packed  and 
unpacked,  put  it  at  the  end  of  the  ix  array. 

ix(len,5)=<  2**31-1> 

print*,'  ***  case  1  ***' 

write  some  packed  fl.  pt.  data  and  read  it  back  in. 
itau  is  .GE.  0,  so  it  is  appended  to  ifile,  itype  is 
.GT.  0,  so  nfwrit  returns  64  bit  data  with  precision 
truncated  to  32  bits. 

itype=  2 

print*,  '  x  before  packing*  ',x(len,5) 

cal l  nfwrit< ifile, 1 , itype, 1 ,n, len,x, itau, cdtg, istat) 

print*,  *  x  returned  from  nfwrit*1,  x<(en,5) 

call  nfread(if i le,1, itype, 1,n, len,x, itau, cdtg, istat) 

print*,  1  x  returned  from  nfread=',  x(len,5) 

print*,'  ***  case  2  ***' 

write  some  packed  integer  data  and  read  it  back  in 
the  data  is  added  to  ifile  starting  at  record  no.  6. 

itype*  1 

print*,  1  ix  before  packing*',  ix(len,5) 

call  nfwrit(ifi le,1,itype,6,n, len, ix, itau, cdtg, istat) 

call  nfread(if i le, 1 , itype,6,n, len, ix, itau, cdtg, istat) 

print*,  '  ix  after  unpacking*' , ix( len, 5) 

print*,'  ***  case  3  ***' 

write  some  packed  fl.  pt.  data  and  read  it  back  in. 
itau  is  .GE.  0,  so  it  is  appended  to  ifile,  itype  is 
.LT.  0,  so  nfwrit  returns  full  64  bit  data,  data  is 
still  written  as  32  bit  IEEE,  however. 

itype*  -2 

print*,  1  x  before  packing*  ',x(len,5) 

cal l  nfwri t( if i le,1, itype, 1,n, len,x, itau, cdtg, istat ) 

print*,  1  x  returned  from  nfwrit*',  x(len,5) 

call  nfread(if i le,1, itype,1,n, len, x, i tau, cdtg, istat) 

print*,  '  x  returned  from  nfread*',  x(len,5) 

print*,'  ***  case  4  ***' 
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c  write  some  packed  integer  data  to  a  file  without  the  appended 
c  counter  and  read  it  back  in 
c 

!  itau=  -1 

!  itype=  1 

!  print*,  1  ix  before  pack ing=' ,  ix(len,5) 

!  call  nfwrit(ifile,1, itype, 1,n, len, ix, itau, cdtg, istat) 

!  call  nfread(if i le, 1 , itype, 1 ,n, len, ix, i tau.cdtg, istat) 
c 

!  print*,  1  ix  after  unpacking=' , ix( len, 5) 
c 

!  print*,'  ***  case  5  ***' 

c 

c  write  some  unpacked  fl.  pt.  data  to  a  file  and  read  it  back  in 
c  note  that  a  different  value  of  positive  itau  generates  a  new 
c  filename  which  can  have  the  different  length  records  the  full 
c  precision  array  requires 
c 

!  itau=  20 

!  i type=  0 

!  print*,  '  xfull  before  write=  ‘,xfull(len,5) 

!  call  nfwrit(ifile,1, itype, 1,n, len, xfull, itau, cdtg, istat) 

!  print*,  1  xfull  returned  f rom  nfwri t=' ,  xfull(len,5) 

!  call  nfread(ifile,1, itype, 1,n,len, xfull, itau, cdtg, istat) 

!  print*,  1  xfull  returned  f rom  nfread=' ,  xfull(len,5) 

c 

c  after  running  this  test  example  look  at  contents  of  your 
c  target  directory  to  verify  filenames  and  sizes, 
c 

!  stop 

!  end 

c 

c  ***END  OF  EXAMPLE******************************************* 
c 

dimension  x(len,nrec) 
dimension  xpack(( len+1 )/2) 
c 

character*8  cdtg 
character*48  fnam,blk24 
character*54  file 
c 

parameter  (nf=20) 
common/fncom/  cfiles(nf) 
conmon/fucom/  fcall 
character*48  cfiles 
logical  fcall, iread 
c 

data  bl k24/ 1 
data  fcal l/. true./ 
if (fcall)  then 
fcal l=. false, 
do  5  i=1,nf 
cfiles(i)=blk24 
5  continue 
endif 
c 

do  10  k=1 ,nf 
kk=  k 

if(cf i les(k).eq.blk24)  cfiles(k)=  fnam 
if (cf i les(k).eq.fnam)  go  to  20 
10  continue 
20  iun=  10+kk 
c 
c 

iread=.true. 

cal l  nfopen(fnam,msg, iread, itau, itype, len, iun.fi le, istat) 


if(istat.eq.2)  then 
print  500,  fnam, istat 

500  format (1x,a54, 1  missing,  istat=',i1) 
return 
endif 
c 

if(itype.eq.O)  then 
c 

c  no  unpacking 
c 

do  30  j=1,nrec 
irec=  istrt+j-1 

read(iun,rec=irec,err=45)  (x(i, j), i=1, (en), itau.cdtg 
30  continue 
e 

else 

c 

c  data  is  packed,  convert  from  ieee  32  bit  to  cray  64  bit 
c 

itp=  abs(itype) 
lenx=  (len+1)/2 
do  40  j=1,nrec 
irec=  istrt+j-1 

read(iun,rec=irec,err=45)  (xpack( i ), i=1 , lenx) , i tau, cdtg 
ierr=  ieg2cray( i tp, len,xpack,0,x(1 , j ) ) 

40  continue 
e 

endif 

c 

if(itau.ge.O)  close(unit=iun) 
c 

istat=  0 

if(msg.eq.O)  return 
c 

print  100, f i le.nrec, istrt, len 

100  formate  lx, a54, 'read:  *,i3,'  records  starting  at  rec=',i4 
:  len=  • , i8> 
c 

if(itype.eq.O)  then 
print  200,  itau.cdtg 

200  formate  data  read  for  tau=',’6.'  :  dtg=  ‘,a8) 
else 

print  300,  i type, itau.cdtg 

300  formate  type', i 2,'  packed  data  read  for  tau=',i6, '  :  dtg=  ',a8) 
endif 
c 

return 

c 

45  print  400,  f i le, iun, irec 

400  formate  bad  read  on  ',354,'  :  unit=,,i3,1  :  record^  ' , i 4 > 
istat=  1 
return 
end 
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subroutine  nfwrit(fnam,msg, itype, istrt, nrec, len,x, itau.cdtg, istat) 
c 

c  subroutine  to  write  NOGAPS  history  files  (but  can  also  be  used  for 
c  general  purpose  10).  file  structure  is  FORTRAN  direct  access,  as  an 
c  option  can  take  CRAY  64  bit  data  and  convert  it  to  32  bit  IEEE  format 
c  (e.g.  SUN  format)  before  writing  the  data,  files  are  opened  internally, 
c  with  up  to  20  files  opened  simulaneously,  using  units  11-30.  depending 
c  on  actual  ncmber  of  files  opened  by  nfread/nfwrit,  users  must  ensure 
c  that  needed  unit  numbers  are  not  used  elsewhere  in  their  program, 
c 

c  formal  parameters 
c 

c  INPUT: 
c 

c  fnam:  character*48  path/filename  to  data  to  be  written,  path  can 
c  relative  or  absolute  but  must  be  pre-existing,  actual  filename  may 

c  have  an  indexing  string  appended  depending  of  value  of  formal  parameter 

c  'itau1.  see  example  below, 

c 

c  msg:  flag  to  control  information  prints  about  file  activity 
c  =0,  no  prints 

c  =1,  print  info  about  read  or  write  action 

c  =2,  print  info  about  file  opening,  reading  and  writing  action 

c  =  anything  else,  same  as  2 

c 

c  itype:  data  type 

c  =  0,  do  not  upack  data,  used  if  file  data  is  already  CRAY  format 

c  =1  (or  -1),  pack  from  CRAY  integer  to  IEEE  integer 

c  =2  (or  -2),  pack  from  CRAY  float  pt.  to  IEEE  float  pt. 

c  if  itype  .GT.  0  then  64  bit  data  words  returned  in  array  x  are  truncated 
c  back  to  32  bit  precision, 

c 

c  istrt:  beginning  direct  access  index  of  records  written  to  fnam 
c 

c  nrec:  number  of  records  (beginning  at  istrt)  written  to  fnam 
c 

c  len:  number  of  data  items  in  array  x  for  fnam  records,  (note:  since 
c  file  format  is  direct  access,  len  must  be  the  same  for  all  records 
c  on  a  given  fnam) 

c 

c  x:  array  containing  data  to  be  written,  should  have  dimension  at  least 
c  (nrec* len) 

c 

c  itau:  counter  index,  if  .GE.  to  zero  it  is  appended  as  a  character*6 
c  string  to  fnam,  yielding  a  new  filename  which  is  the  actual 
c  one  opened,  itau  is  also  written  out  as  part  of  data  record, 
c  see  example  below, 

c 

c  cdtg:  data- time-group  character*8  item  of  data  written  as  part  of  data 
c  record 

c 

c  OUTPUT: 
c 

c  istat:  returned  status  code, 
c  =0,  no  problems 

c  =1,  bad  write,  data  not  output 

c 

c  EXAMPLE: 

c 

£  ******************************************************************** 
c 

!  program  test 

c 

c  program  to  demonstrate  use  of  subroutines  nfread  and  nf writs 
c  remove  Ms  from  colunn  1  to  get  compilable  code 
c 

!  parameter  ( len=  10001) 
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c 

c 

c 

! 

j 

! 

! 

! 

j 

! 

c 

c 

c 

c 

! 

c 

! 

c 

c 

c 

c 

c 

c 

! 

! 

i 


c 

I 

c 

c 

c 

c 

I 


! 

c 

I 

c 

! 

c 

c 

c 

c 

c 

c 


dimension  x(ten,5) 
dimension  xfull(len,5) 
dimension  ix(len,5) 
character*48  ifile 
character*8  cdtg 

change  the  following  path/filename  to  suit  your  interests 

ifi  le='/a/yourname/datapath/dunmy ' 

cdtg='92030412' 

n=5 

i tau=  10 

generate  some  test  data 

do  1  j=1,n 
do  1  i=1,len 
x(i, j)=  (i+j) 
x(i, j)=-x(i, j >**3.888 
xfullfi, j)=  x( i , j ) 
ix( i , j )=  ( i+j )**2 
1  continue 

this  is  the  largest  integer  preserved  when  packed  and 
unpacked,  put  it  at  the  end  of  the  ix  array. 

ix(len,5)=(  2**31-1) 

print*,1  ***  case  1  ***' 

write  some  packed  fl.  pt.  data  and  read  it  back  in. 
itau  is  .GE.  0,  so  it  is  appended  to  ifile,  itype  is 
•GT.  0,  so  nfwrit  returns  64  bit  data  with  precision 
truncated  to  32  bits. 

i type=  2 

print*,  '  x  before  packings  1 ,x(len,5) 

cal l  nfwri t( if i le, 1 , i type, 1 ,n, len,x, itau, cdtg, istat) 

print*,  1  x  returned  from  nfwri t= 1 ,  x(len,5) 

cal l  nf read(if i le, 1 , itype, 1 ,n, len,x, i tau.cdtg, istat) 

print*,  1  x  returned  from  nfreads',  x(len,5) 

print*,1  ***  case  2  ***' 

write  some  packed  integer  data  and  read  it  back  in 
the  data  is  added  to  ifile  starting  at  record  no.  6. 

i type=  1 

print*,  1  ix  before  packings',  ix(len,5) 

call  nfwri t(if i le, 1, i type, 6, n, len, ix, i tau.cdtg, istat) 

call  nfread< if ile, 1, i type, 6, n, len, ix, i tau.cdtg, istat) 

print*,  1  ix  after  unpackings 1 , ix( len, 5) 

print*,'  ***  case  3  ***' 

write  some  packed  fl.  pt.  data  and  read  it  back  in. 
itau  is  .GE.  0,  so  it  is  appended  to  ifile,  itype  is 
.LT.  0,  so  nfwrit  returns  full  64  bit  data,  data  is 
still  written  as  32  bit  IEEE,  however. 

itype=  -2 

print*,  '  x  before  packings  >,x(len,5) 

cal l  nfwri t( if i le, 1 , itype, 1 ,n, len,x, i tau.cdtg, istat) 

print*,  '  x  returned  from  nfwrit=',  x(len,5) 

cal l  nf read( if i le, 1, i type, 1,n, len,x, i tau.cdtg, istat) 

print*,  1  x  returned  from  nfread=',  x(len,5) 


A-l  5 


c 

!  print*,1  ***  case  4  ***' 
c 

c  write  some  packed  integer  data  to  a  file  without  the  appended 
c  counter  and  read  it  back  in 
c 

I  i tau=  - 1 

I  i type=  1 

•  print*,  '  ix  before  packings' ,  ix(len,5) 

•  call  nfwrit(ifile,1,itype,1,n,len,ix,itau,cdtg,istat) 

!  call  nfread(ifile,1,itype,1,n,len,ix,itau,cdtg,istat) 

c 

I  print*,  '  ix  after  unpackings' , ix( len, 5) 
c 

!  print*,'  ***  case  5  ***' 

c 

c  write  some  unpacked  fl.  pt.  data  to  a  file  and  read  it  back  in. 
c  note  that  a  different  value  of  positive  itau  generates  a  new 
c  filename  which  can  have  the  different  length  records  the  full 
c  precision  array  requires 
c 

!  i tau=  20 

!  i type=  o 

!  print*,  1  xfull  before  writes  ' ,xful l( len, 5) 

!  call  nfwrit(ifile,1,itype,1,n, len, xfull, itau, cdtg.istat) 

!  print*,  '  xfull  returned  from  nfwrit=‘,  xfull(len,S) 

!  call  nfread(ifile,1,itype,1,n, len, xfull, itau, cdtg.istat) 

!  print*,  •  xfull  returned  from  nfread=' ,  xfull(len,S) 
c 

c  after  running  this  test  example  look  at  contents  of  your 
c  target  directory  to  verify  filenames  and  sizes, 
c 

!  stop 

!  end 

c 

c  ***enp  OF  EXAMPLE******************************************* 
c 

dimension  x(len,nrec) 
dimension  xpack(( len+1 )/2) 
c 

character*8  cdtg 
character*48  fnam,blk24 
character*54  file 
c 

parameter  (nf=20) 
common/fncom/  cfiles(nf) 
common/ fucom/  fcall 
character*48  cfiles 
logical  fcall, iread 
c 

data  blk24/'  '/ 

data  fcal 1/ . true./ 
if(fcall)  then 
fcal l=. false, 
do  5  i=1,nf 
cf i les( i )=blk24 
5  continue 
endi  f 
c 

do  10  k=1,nf 
kk=  k 

if(cf i les(k).eq.blk24)  cfiles(k)=  fnam 
if(cfiles(k).eq.f nam)  go  to  20 
10  continue 

20  iuns  10+kk 
c 

itp=  abs(itype) 
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iread=. false. 

call  nfopen(fnam,msg, iread, i tau, i  tp, len, iun.fi le, istat) 
c 

if(itype.eq.O)  then 
c 

c  write  data  as  is,  no  packing 
c 

do  30  j=1,nrec 
irec=  istrt+j-1 

write(iun,rec=irec,err=45)  <x(i,j),i=1, len), itau.cdtg 
30  continue 
c 

else 

c 

c  pack  data,  cray  64  bit  to  ieee  32  bit 
c 

lenc=  (len+1)/2 
itp=  abs(itype) 
do  40  j=1,nrec 
irec=  istrt+j-1 

if(itype.eq.-2)  call  undf lo(x(1 , j ), ten) 
ierr=  cray2ieg(itp, len,xpack,0,x(1 , j ) ) 
write(iun,rec=i rec,err=45)  (xpack( i ), i=1 , lenc), itau.cdtg 
if(itype.gt.l)  ierr=  ieg2cray( i tp, len, xpack,0,x( 1 , j ) ) 

40  continue 
c 

endi  f 
c 

if(itau.ge.O)  close(unit=  iun) 
c 

istat=  0 

if(msg.eq.O)  return 
c 

print  100, f ile.nrec, istrt, len 

100  formate  lx, a54, 'write:  * , i3, ’  records  starting  at  rec-', i4 
*,'  :  len=  ',i8) 
c 

if(itype.eq.O)  then 
print  200,  itau.cdtg 

200  formate  data  written  for  tau=',i6,'  :  dtg=  ',a8) 
else 

print  300,  itype, itau.cdtg 

300  formate  type',i2,1  packed  data  written  for  tau=',i6,'  :  dtg=’,a8) 
endif 
c 

return 

c 

45  print  400,  f i le, iun, i rec 

400  formate  bad  write  on  ',a54,'  :  unit=',i3,1  :  record=  * ,  i 4 ) 
istat=  1 
return 
end 
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subroutine  noday(cidtg,my,mhy) 
c 

c  subroutine  to  build  date-time-group  from  hour  of  year 
c 

c  input 

c  my  =  last  2  digits  of  year 

c  mhy=  hour  of  the  year 

c 

c  output 
c 

c  return  ascii  dtg  (yymmddhh)  in  idtg 
c 

c  ****••*••**••*•*••**••*•**•**•**•***••********•***•**•**** 
c 

character*8  cidtg 
dimension  month(12,2) 

data  month/0, 31, 59, 90, 120, 151 ,181 ,212, 243,273,304, 334, 
S  0,31,60,91,121,152,182,213,244,274,305,335/ 

c 

j  =  1 

mdy  =  mhy/24+1 
if(mod(my,4).eq.0)  j  =  2 
kk-0 

do  2  i  =  2,12 
kk=kk+1 

if (mdy. le.month( i , j))  go  to  3 

2  continue 
kk=kk+1 

3  mm  =  kk 

mh  =  mod(mhy,24) 
md  =  mdy-month(mm, j) 
write(cidtg,1)  my,mm,md,mh 
1  format(4i2) 
do  5  i  =  1,8 

if (cidtg(i : i ).eq. 1  ■)  cidtg< i : i )='0' 

5  continue 

c 

return 

end 
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subroutine  sortmK jtrun, mlmax, msort, Isort, mlsort) 
c 

c  subroutine  to  generate  sorting  arrays  for  locating  elements 
c  in  NOGAPS  spherical  harmonics  arrays 
c 

c  i nput : 

c 

c  jtrun:  truncation  limit:  =48  for  T47  model 
c  mlmax:  ( jtrun+1)*jtrun/2 
c 

c  output: 

c 

c  msort:  zonal  waveno(m)=  msort(ml),  mt=1, mlmax 
c  Isort:  total  waveno(l)=  Isort(ml),  ml=1(mlmax 
c  mlsort:  triangular  truncation  spherical  harmonic 

c  index(ml)=  mlsort(m,l) 

c 

c  **************************************************** 

c 

dimension  msort (mlmax), Isort (mlmax) , ml  sort ( jtrun, j t run) 
c 

mlx=  ( jtrun/2)*((jtrun+1)/2) 
ml=  0 

do  1  k=1, jtrun-1,2 

do  1  m=1,jtrun-k 

ml=  ml  +  1 
mlp=  ml+mlx 
mlsort(m,m+k)=  mlp 
mlsort(m,m+k-1 )=  ml 
msort(ml)=  m 
lsort(ml)=  m+k-1 
msort(mlp)=  m 
lsort(mlp)=  m+k 

1  continue 
c 

ml=  mlp 

do  2  m=2,jtrun,2 
ml=  ml+1 

mlsort(m, jtrun)=  ml 
msort(ml)=  m 
lsort(ml)=  jtrun 

2  continue 
return 
end 


A-  1  9 


subroutine  tothrs< idtg, ichour) 
c 

c  subroutine  converts  dtg(idtg)  to  hour  of  century, 
c 

c  input: 
c 

c  idtg:  date-time-group  (YYMMDDHH) 
c 

c  output : 

c 

c  hours  since  midnight  1  Jan  1900 
c 

C  ********************************************************** 

c 

characters  idtg.icdtg 
read( idtg, 1 (4i2) 1 )iyy, imn, idd, ihh 
icdtg=idtg 
iehour=0 
iyear=0 
c 

c  compare  for  current  year  .jump  to  final  calculations, 
c 

100  if(iyear  .eq.  iyy)  go  to  500 
jday=365 
c 

c  check  to  year  at  00  .  1900  uas  not  a  leep  year  so  leep  test 
c  isn't  executed. 

c  compare  current  year  to  override  leap  year  check, 
c  allows  to  subroutine  to  execute  till  the  year  2040. 
c  ****  please  note  changing  below  line  will  allow  subroutine 
c  to  execute  longer  of  shorter. depending  of  interger  given, 

c  in  the  year  2100  deleting  the  following  two  lines  will  allow 

c  this  subroutine  to  run  till  the  year  2400. 

c 

if(iyy  .It.  40)  goto  300 
c 

c  check  for  leap  year, if  yes,  add  1  day  to  jdaysfdays  in  the  year) 
c 

300  if (mod( iyear,4)  .eq.  0)  jday=jday+1 

200  ihour=jday*24 
iyear  =  iyear+1 
i chour= i chour+ i hour 
go  to  100 
c 

c  calculate  current  years  hour, 
c 

500  call  daynum( icdtg, iyy, imm, idd, ih, idanum, ihh) 
ichour=ichour+ihh 
return 
end 
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subroutine  tranrs( istrt, l l ,poly,w, r,s) 
c 

c  subroutine  to  transform  T 47  NOGAPS  gaussian  grid  fields 
c  to  spherical  harmonic  arrays 
c 

c  i nput : 
c 

c  istrt:  quadrature  initialization  flag 
c  eq  0,  initialize  array  's'  from  zero 

c  ne  0,  add  contributions  to  pre-existing  values  in  's' 

c  ll:  number  of  levels( layers)  of  input/ouput  arrays 
c  poly:  associated  legendre  polynomials  (from  subroutine  Igndr) 
c  w:  gaussian  quadrature  weights  (from  subroutine  gausl3) 
c  r:  input  array  of  gaussian  grid  fields  to  transform 
c 

c  output : 

c 

c  s:  output  array  of  spherical  harmonic  coeff  fields 
c 


c 

parameter  (im=144,  lm=18,  lpout=6) 
parameter  (jm=  im/Z,  im2-  im/2,  jm2=  jm/2) 

parameter  (jtrun=  2*{ ( 1+( im- 1 )/3)/2),  mlmax=  ( jtrun+1 )*jtrun/2) 
parameter  (imlm=im*lm,  imjm=im*jm,  idim2=  mlmax*2) 
c 

parameter  (mlx=  (jtrun/2)*(( jtrun+1 )/2)) 
parameter  (junp=  im+3) 
c 

dimension  w( jm) ,poly( idim2, jm2) ,s(mlmax,2, ll),r<im,jm,ll) 
dimension  cc( jump, jm) ,uork( imjm,2) 
c 

cornnon/  fftcom/  trigs( im,2) , i fax( 19) 
c 

jstrt=  1 

if(istrt.eq.O)  jstrt=  2 
c 

do  30  k=1 . 1 1 
do  23  j=1,jm 
do  23  i  =  1 , im 
cc(i,j)=  r(i, j,k) 

23  continue 
c 

cal  l  rfftmlt(cc, work, trigs,  ifax,  1 ,  jump,  im,  jm,  -1 ) 
c 

c  if  istrt  .eq.  zero,  the  quadrature  integral  is  initialized  from  zero, 
c  otherwise  the  sum  is  added  to  initial  's'  array  passed  with  call, 
c 

i f ( i strt .eq.O)  then 
c 

m1=  0 

do  62  l=jtrun-1,1,-2 
CDIR3  IVDEP 

do  63  m  =  1,  l 
mm=  2*m-1 
mp=  mm+1 
ml=  m+ml 
mk=  ml+mlx 

s(ml,1,k)  «  w(1)*poly(ml,1)*(cc(mm,1)+cc(irm,  jm)) 
s(ml,2,k)  =  w(1 )*poly(ml , 1 )*(cc(mp, 1 )+cc(mp, jm)) 
s(mk,1,k)  =  w(1)*poly(mk,1)*(cc(mm,1)-cc(nn, jm)) 
s(mk,2,k)  =  w(1)*poly(mk,1)*(cc(mp,1)-cc(mp, jm)) 

63  cont i nue 

m1=  ml+l 
62  continue 
c 

ml=  mlx*2 
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COIRS  IVDEP 

do  64  m=2,jtrun,2 
mt=ml+1 
«nt=  2*m-1 
mp=  mm+1 

s(ml , 1 ,  k)=  w(1 )*poly(ml , 1 )*(cc(mm, 1 )+cc(mm, jm) ) 
s<ml ,2, k)=  w( 1 )*poly(ml , 1 )*(cc(mp, 1 )+cc(mp, jm) ) 
64  continue 
endif 


c 


do  70  j=jstrt,jm2 
jj=  jm-j+1 
m1=  0 


COIRS 


73 


do  72  l=jtrun-1,1,-2 


IVDEP 

do  73  m  =  1, 
mm=  2*m-1 
mp=  mm+1 
ml=  nt+ml 
mk=  ml+mlx 
s(ml,1,k) 
s(ml,2,k) 
s(mk,1,k) 
s(mk,2,k) 
cont i nue 


s(mt , 1 ,k)+w( j )*poly(ml , j )*(cc(nm, j >+cc(mn, j  j  > ) 
s(ml ,2,k)+u( j )*poly(ml , j )*(cc(mp, j )+cc(mp, j  j )) 
s(mk,  1  ,k)+n(  j)*poly(mk,  j  )*(cc(mm,  j  )-ccCmm,  j  j )) 
s<mk,2,k)+w( j)*poly(mk, j)*(cc(mp, j)-cc(mp, j j)) 


m1=  ml* l 


72  continue 


c 

ml=  mlx*2 


COIRS  IVDEP 

do  65  m=2, jtrun,2 

ml=ml+1 

nm=  2*m-1 


mp=  mm+1 

s(ml , 1 , k)=  sCml , 1 ,k)+w( j)*poly(ml , j )*(cc(im,  j )+cc(mm,  j  j )) 
s(ml ,2,k)=  s(ml,2,k)+w( j )*poly(ml , j )*(cc(mp, j )+cc(mp, j  j )) 
65  continue 
70  continue 
30  continue 


return 

end 
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subroutine  transr(ll,poly,s,r) 


subroutine  to  transform  a  T47  NOGAPS  spherical  harmonic  coefficient 
arrays  to  equivalent  gaussian  grid  gridpoint  arrays 

input: 

l l :  number  of  levelsf layers)  in  input/output  arrays 

poly:  associated  legendre  polynomials  (from  subroutine  Igndr) 

s:  spherical  harmonic  arrays 

output: 

r:  gaussian  grid  arrays 

************************************************************** 

parameter  <im=144,  lm=18,  lpout=6) 
parameter  (jm=  im/2,  im2=  im/2,  jm2=  jm/2) 

parameter  (jtrun=  2*((1+( im- 1 )/3)/2) ,  mtmax=  ( jtrun+1 )*jtrun/2) 
parameter  (imlm=im*lm,  imjm=im*jm,  idim2=  mlmax*2) 

parameter  (mlx=  ( j trun/2)*( ( j trun+1 )/2)) 
parameter  (jump=  im+3) 

dimension  poly( idim2, jm2) ,s(mlmax,2, l l ), r( im, jm, l l) 
dimension  cc( jump, jm),uork( imjm,2) 

conmon/  fftcom/  trigs( im,2), i fax(19) 


do  20  k=1,tl 
do  55  m=1,jump*jm2 
cc(m,1)=  0.0 
cc(m, jm2+1)=  0.0 
55  continue 


do  5  j=1,jm2 
jj=  jm*1  -  j 
ml=  2*mlx 
CD  IRS  IVDEP 

do  3  m=2, j trun, 2 
ml=  ml+1 
mm=  2*m-1 
mp=  mm*1 

cc(mm,j)=  polyfml, j)*s(ml,1,k) 
cc(mp,j)=  poly(ml , j )*s(ml ,2, k) 
cc(mm,jj)=  cc(mm, j ) 
cc<mp, j  j )=  cc(mp, j ) 

3  continue 


c 

ml-  0 

do  5  l=jtrun-1,1,-2 
CD  IRS  IVDEP 

do  6  m=1, l 
mm=  2*m- 1 
mp=  iun*1 
ml=  m*m1 
mk=  ml+mlx 

cc(mn,  j  )=  ccCnrm,  j  )+poly(ml ,  j  )*s(ml  ,1  ,k)+poly(mk,  j  )*s(mk,  1  ,k) 
cc(rrm,  jj)*cc(mni,  jj)+poly(mt,  j)*s(ml,1,k)-poly(mk,  j)*s(mk,1,k) 
cc(mp, j)=  cc(mp, j )+poly(ml ,j)*s(ml,2,k)+poly<mk,j)*s(mk,2,k) 
cc(mp, j j)=cc(mp, jj)+poly(ml, j )*s(ml ,2, k)-poly<mk, j)*s(mk,2,k) 
6  continue 
m1=  ml  +  l 
5  continue 


c 

cal l  rfftmltfcc, work, trigs, i fax,1 , jump, im, jm, 1) 
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do  22  j=1 
do  22  i=1 
r(i, j,k)= 
22  continue 
20  continue 

return 

end 


,  jm 
,  im 

cc(i, j) 


A- 24 


subroutine  tranuv( l l ,radsq,oeos,eps4,cim,poly,dpoly 
*,  vor,div,ut,vt) 
c 

c  subroutine  to  transform  NOGAPS  T47  spherical  harmonic  arrays  of 
c  vorticity  and  divergence  into  gaussian  grid  u  and  v  components 
c 

c  input: 
c 

c  ll:  number  of  levels( layers)  of  input  and  output  fields 
c  radsq:  (radius  of  earth)**2  =  (6371000.0**2) 
c  ocos:  1 .0/(cos( l at i tude)**2 

c  eps4:  spherical  harmonic  laptacian  operators  l*( 1+1 )/radsq,  where 
c  1=  total  waveno  ( 1=0, jtrun-1 )  (from  subroutine  cons) 

c  cim:  zonal  wavenunber  operator  (from  subroutine  cons) 
c  poly:  associated  legendre  polynomials  (from  subroutine  Igndr) 
c  dpoly:  d(poly)/d(sin( lati tude) )  (from  subroutine  Igndr) 
c  vor:  spherical  harmonic  vorticity 
c  div:  spherical  harmonic  divergence 
c 

c  output : 
c 

c  ut:  u-component  (weighted  by  cos(latitude)/(earth  radius)) 
c  vt:  v-component  (weighted  by  cos( lat i tude)/(earth  radius)) 
c 

c  ************************************************************** 

c 

parameter  (im=144,  lm=18,  lpout=6) 
parameter  (jm=  im/2,  im2=  im/2,  jm2=  jm/2) 

parameter  (jtrun=  2*(( 1+( im-1 )/3)/2) ,  mlmax=  ( jtrun+1 )+ j trun/2) 
parameter  (imlm=im*lm,  imjm=im+jm,  idim2=  mlmax*2) 
c 

parameter  (mlx=  ( jtrun/2)*( ( jtrun+1 )/2) ) 
parameter  (jumps  im+3) 
c 

dimension  ocos( jm),eps4(mlmax),cim(m!max),poly( idijn2, jm2) 

*,  dpoly(idim2, jm2),vor(mlmax,2, Im) 

*,  div(mlmax,2, Im) ,ut( im, jm, Im) , vt < im, jm, Im) 
c 

dimension  cu( jump, jm),cv( jump, jm),work( imjm,2) 

*,  cfac(mlmax),dfac(mlmax) 
save  cfac.dfac 
c 

common/  fftcom/  trigs( im,2), ifax(19) 
c 

logical  first 
data  first/. true./ 
c 

if(first)  then 
first=. false. 

CDIRa  IVDEP 

do  10  mt=2,mlmax 

dfac(ml)=  1.0/(radsa*eps4(ml)) 

cfac(ml)=  cim(ml)*dtac(ml ) 

10  continue 
dfac(1)=  0.0 
cfac(1)=  0.0 
endi  f 
c 

do  15  k=1,ll 
do  55  m=1,jump*jm2 
cu(m, 1)=  0.0 
cu(m, jm2+1)=  0.0 
cv(m,1)=  0.0 
cv(m,jm2+1)=  0.0 
55  continue 
c 

do  20  j=1 , jm2 
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c 


rcos=  1.0/ocos(j) 

jj*  jm+1 - j 


ml=  2*mlx 
CDIR3  1VDEP 

do  50  m=2,jtrun,2 
ml=  ml+1 
2*m-1 
mp=  wn+1 

cu(rm>,  j  )=  ♦cf  ac(mt  )*poly(ml ,  j  )*di  v(mi ,2, k) 

*  +df ac(ml )*dpoly(mt , j )*rcos*vor(mt , 1 , k) 
c 

cuCmm, j j )=  ♦cfac(ml )*poly(ml , j)*div(ml,2,k) 

*  -df ac(mt )*dpot y(mt , j )*rcos*vor(ml , 1 ,  k) 


c 

cu(mp, j )=  -cf acCml )*poty(ml , j )*di v(ml , 1 , k) 

*  +df ac(ml )*dpoly(ml , j )*rcos*vorCmt ,2, k) 


c 


c 


cu(rap, jj)= 

* 


cvCmm, j )= 


c 


c 


c 


c 


cvCmm, j j )= 

* 

cv(mp, j )= 

* 

cv(mp, jj)= 

* 

50  continue 


-cfac(ml  )*pot y(ml , j )*di v(ml ,  1 , k) 

-df ac(ml )*dpoly(ml , j )*reos*vor(mt ,2, k) 

♦cfacCml )*poly(ml , j)*vor(ml ,2, k) 
-dfac(ml )*dpoty(ml , j )*rcos*di v(mt ,  1 ,  k) 

♦cfacCml )*poly(ml , j )*vor(ml ,2,k) 
♦dfacCml )*dpoly(ml , j )*rcos*di v<mt , 1 ,k) 

-cfacCml )*poly(ml ,j )*vor(ml , 1 ,k) 
-dfacCml )*dpoly(ml , j )*rcos*di vCml ,  2,k) 

-cfacCml  )*polyCml ,  j  )*vor(ml ,  1 ,  k) 
♦dfacCml )*dpolyCml, j)*rcos*div(mt,2,k) 


m1=  0 

do  30  1= jtrun-1 , 1,-2 
COIRS  IVDEP 

do  40  m=1 , l 
irm=  2*m-1 
mp=  rrrn+1 
ml  =  m*m1 
mk=  ml+mlx 


c 

cuCmm,  j  )=  cuCmm,  j  )+cfac(ml  )*polyCml ,  j  )*di  v(ml  ,2,  k) 

*  +cfacCmk)*polyimk, j )*di vCmk,2,k) 

*  +dfac(ml )*dpolyCml , j )*rcos*vorCml , 1 , k) 

*  +dfacCmk)*dpolyCmk, j )*rcos*vor(mk, 1  ,k) 
c 

cuCmm, j  j )=  cu(mm, j  j )+cfac(ml )*polyCml , j )*di vCml ,2, k) 

*  -cfacCmk)*polyCmk, j)*divCmk,2,k) 

*  -dfac(ml )*dpolyCml , j )*rcos*vor(ml , 1 , k) 

*  +dfacCmk)*dpolyCmk, j )*rcos*vor(mk, 1 , k) 
c 

cu(mp, j)=  cuCmp. j)-cfac(ml )*polyCml, j )*di vCml , 1 ,k) 

*  -cf ac(mk)*poly(mk, j )*di vCmk,1 ,k) 

*  +dfac(ml )*dpolyCml, j)*rcos*vor(rol ,2,k) 

*  ♦dfacCmk)*dpolyCmk, j )*rcos*vor(mk,2,k) 
c 

cuCmp, jj)=  cuCmp, j j ) -cfacCml )*polyCml , j)*divCml,1,k) 

*  +cfacCmk)*polyCmk, j )*div(mk, 1 ,k> 

*  -dfacCml )*dpolyCml , j )*rcos*vorCml ,2, k) 

*  ♦dfacCmk)*dpolyCmk, j )*rcos*vor<mk,2,k) 
c 

cvCmm,  j  )=  cvCmm,  j  )+cf acCml  j*polyCml ,  j  )*vorCml  ,.-,k) 

*  ♦cfacCmk)*polyCmk,  j )*vor(mk,2,k) 

*  -dfacCml )*dpotyCml , j )*rcos*di v<ml , 1 , k  > 
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*  -dfac(mk)*dpoly(mk, j )*rcos*di v(mk, 1 , k) 
c 

cv(n*n,  j  j  )=  cv(mm,  j  j  )+cfac(ml  )*poly(ml ,  j  )*vor(ml ,  2,k) 

*  -cfac<mk)*poly(mk, j)*vor(mk,2,k) 

*  +dfac(ml )*dpoly(ml , j )*rcos*di v(ml , 1 , k) 

»  -dfac<mk)*dpoly(mk, j )*rcos*div(mk, 1 ,k) 

c 

cv(mp, j )=  cv(mp, j )-cfac(ml )*poty(ml , j )*vor(ml , 1 , k) 

*  -cfac(mk)*poty(mk, j )*vor(mk, 1 ,k) 

*  -dfac(ml)*dpoly(mt , j)*rcos*di v<ml ,2,k) 

*  -dfac<mk)*dpoly(mk( j )*rcos*di v(mk,2,k) 
c 

cv(mp, j  j )=  cv(mp, j  j )-cf ac(mt )*poly(ml , j )*vor(ml , 1 , k> 

*  +cfac(mk)*poty(mk( j)*vor(mk, 1,k) 

*  +dfac(ml )*dpoly(ml , j )*rcos*di v(ml , 2, k) 

*  -dfac(mk)*dpoly(mk, j )*reos*di v(mk, 2, k> 
c 

40  continue 
c 

m1=  ml+l 
30  continue 
20  continue 
c 

call  rfftml t(cu,uork, trigs, i fax, 1 , jump, im, jm, 1) 
cal l  rf ftmlt(cv, work, trigs, ifax, 1 , jump, im, jm, 1 ) 
c 

do  22  j=1,jm 
CDIR3  IVDEP 

do  22  i=1 , im 
ut(i, j,k)=  cu( i , j ) 
vt( i , j ,k)=  cv< i , j ) 

22  continue 
c 

15  continue 
return 
end 


A-27 


subroutine  undflo(x,m) 
c 

c  subroutine  to  correct  for  illegal  floating  point  values 
c  generated  when  CRAY  software  converted  underflowed  IEEE 
c  32  bit  values  back  to  64  bits 
c 

c  input 

c 

c  x:  input  64  bit  array  which  may  contain  illegal  fp  values 
c  m:  number  of  elements  in  x 
c 

c  output 

c 

c  x:  corrected  values  of  x 
c 

c  ********************************************************** 
c 

dimension  x(m) 
c 

do  10  i=1,m 

if(abs(x(i)).lt.1.0e-40)  x(i)=  0.0 
10  continue 
return 
end 
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c 


c 

c 


c 


c 


c 


c 


c 


c 

c 

c 

c 


c 


c 


program  makeave 

parameter  (im=144,  jm=im/2,  lm=  18,  lpout=6) 

parameter  (jtrun=  2*((1+(im-1)/3)/2),  mlmax=  ( jtrun+1 )*jtrun/2) 

parameter  (idim2=  2*mtmax) 

parameter  (nspec=4*lm+1 ) 

parameter  (numave=  40+21*lpout+2+14*lpout+7) 

parameter  (imjm=im* jm) 


dimension  polaU<mlmax,2,  jm/2) 

dimension  poty( idim2, jm/2), dpoty( idim2, jm/2) , we i ght( jm> 
dimension  eps4(mlmax),cim(mlmax) ,sinl ( jm) ,cosl( jm) ,msort(mlmax) 

*,  lsort(mlmax),mlsort( jtrun, jtrun),ocos( jm),onocos( im, jm) 

common/  fftcom/  trigsC im,2), ifax(19) 

equivalence  (polal 1(1 , 1 , 1 ),poty) 
equivalence  (polal 1(1 ,2, 1 ) ,dpoly) 

dimension  x( im* jm) ,pt( im, jm),avsum(im*jm),xt(im*jm,36) 

*,  xl1(idim2,nspec),xl2(idim2,nspec),xl3( idim2,nspec) 

*,  xl4(idim2,nspec),slpx(imjm),crf iii(imjm),dswi ii(imjm) 

*,  olri i i(imjm),dlpl(imjm),dtpl(imjm) 

*,  phi (imjm,  Im),  tt(  imjm,  lrrn-1 ), rvor( imjm,  tm),  rdi  v( imjm,  Im) 

*,  rstrm( imjm, lm),vpot( imjm, tm) 

*,  sht(imjm, lm),ut(imjm, lm),vt(imjm, Im) ,sgeo( imjm) 

*,  plog( imjm, lm»1 ),pl l( imjm, lm),pk( imjm, Im) ,pk2( imjm, Im) 

common/pvar/  phips( imjm, Ipout), temps( imjm, Ipout) 

*,  sphum( imjm, lpout),psi ( imjm, Ipout ), chi ( imjm, Ipout) 

*,  upres( imjm, Ipout), vpres( imjm, lpout),slp(imjm) 

* ,  phipsv( imjm, Ipout), tempsv( imjm, Ipout) 

*,  sphumv(imjm, lpout),psiv(imjm, lpout),chiv( imjm, Ipout) 

*,  upresv(imjm, lpout),vpresv( imjm, Ipout).slpv(imjm) 

common/spec/vornow(mlmax,2, lm),divnow(mlmax,2, Im) 

*,  temnow(mlmax,2, lm),shnou(mlmax,2, lm),plnow(mlmax,2) 

equivalence  (xll.vornow) 

common/c3aves/  ave2d( imjm, 36) ,ave3d( imjm, tpout*21) 

common/trp/  pdi f f ( imjm) , tsave( imjm) , tmppl ( imjm) 

common/cross/  tempcs(imjm, Im) ,shtcs( imjm, Im) , rhcs( imjm, Im) 

*,  cldscs( imjm, Im) ,uucs( imjm, Im) , vvcs( imjm, Im) ,strmcs( imjm, Im) 

*,  cstt( jm, lm),cssht( jm, lm),csrh( jm, Im) 

*,  csclds(  jm,  lm),csuu(  jm,  Im) ,csw(  jm,  Im)  ,csstrm(  jm,  Im) 

*,  pcs(lm),ptzm(jm) 

dimension  glob(288*145),xin(288*145),yin(288*145) 

*,  alat(jm*2) 
dimension  vuk(288*145*10) 
dimension  gf l ip( 144*73) 

parameter  (lccn=  350+tm*(12+5*lm)) 

******************************************************************** 


dimension  ccn(lccn) 
equivalence  (c,ccn) 
character*8  idtg 

common/cnstnt/  c(300), i tau, i thist , idtg(3),sarray(lm, Im) 
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c 

c 

c 


c 


c 


c 

c 

c 

c 

c 

c 


c 


c 


*,  carray( lm, lm),evecin( lm, lm),evectr( lm, lm), tmcor( lm, Im) 

*,  eigval(lm),spalm(lm),pmcor(lm),pbcor(lm),tmean(lm) 

*,  asig( lm*1  ),bsig( lm+1  ),dasig( lm),dbsig( lm) 

*,  sige(lm+1),dsig(lm),pkfac(lm) 

*,  ptmean,ccpad(100) 

******************************************************************* 


equivalence 
1  (c(1),jmm1 >, 

2,  <c(5), iday), 

3,  (c(7),itauo), 

4,  (c(10), jsplit), 

5,  (c(13),taue), 

6,  (c(16),cosd), 

7,  (c(19),eccld), 

8,  (c(22),prevap), 

9,  (c(25), omega), 

1,  (c(28),pi2), 

2,  (c(31),cp), 

3,  (c(41),f luxcl), 

4,  <c(44),tauh), 

5,  (c(47),mombdg), 

6,  (c(50),p00), 

7,  (c(53),hybrd), 

8,  (c(56),dgtime), 

equivalence 
1  (c(68) , taud) , 

2,  (c(70),pcon), 

3,  (c(73),skew), 

4,  (c(76), i thfg), 

5,  (c(79),tseres), 

6,  (c(93),daypyr), 

7,  <c(96), iqnx), 

8,  (c(99),eccn), 

9,  <c(102),critl), 

1,  <c(105),pcp), 

2,  (c(109), ipfcst), 

3,  (c(112), ipoutp), 

4,  (c(115),icnmax), 

5,  (c(118),nmodev), 

6,  (c(121),uvmax), 
8,  (c(135),prnt), 

equivalence 
1  (c(138),novar), 

2,  (c(141), nozone), 

3,  (c(144) , incrup), 

4,  < c< 1 5 1 ),nsdedy), 

5, 


(c(3), jmm2), 

(c(6) , tofday) 

(c<8), julian), 

(c(11),taui), 

(c(14),tice), 

(c(17),ddraft), 

(c(20),rbear), 

<c(23),radsq) 

(c(26),tauave), 

(c(29),rgas), 

(c(32),hltm), 

(c(42), isunfx), 

(c(45),krnit), 

(c(48),evaprh), 

(c(51),dta), 

(c(54),ksgeo), 

(c(57),ecadv) 


(c(69),dt) 

(c(71 ) .pconkt ) , 
(c(74),grav2), 
(c(77), itaufg), 
(c(80),prftop), 
(c (94), day), 
(c(97),perhl ), 
(c(100),rad), 
(c(103),p608), 
(c(106),f rqdib), 
(c(1 10), ipsigp), 
(c(113),tau), 
(c(116),frqlim), 
(c(119), dinit), 
(c(123),facrad), 
(c( 136) .capapl ) , 


(c(4), Increc) 

(c(9) , jskip) 
(c(12),dtrad) 
(c(15),forHrd) 
(c(18), jtrunp) 
(c(21),prntij) 

(c(27),pi) 
(c(30),capa) 
(c(33),f im) 
(c(43),ptop) 
(c(46), Is) 
(c(49),vis) 
(c(52),ops) 
(c(55),pok) 


(c(72) , Isl ) 
(c(75),delmf ) 
<c(78), Isx) 
(c(81),dtgpcm) 
(c(95),rotper) 
(c(98),decmax) 
(c(101 ),grav) 
(c(104), Imid) 
(c(107),cdtgoi ) 
<c(111), ipinit) 
(c(114),drgdry) 
(c( 1 17) .beta) 
(c(120), taudb) 
(c(134),gamfac) 
(c( 137) ,cuprh) 


(c(139),cdet),  (c(140) ,dcapa) 

(c(142),sstclm),  (c(143) ,updayc) 

<c< 145 ) , incwnd),  (c( 146) ,pcmdi ) 
(c(152),kmonth),  (c(153),mnthdy) 

(c(157),rsdist),  <c(158),sind) 


:************ 


♦♦♦★♦♦♦♦★★♦★★★★★★★★♦★★★★★♦★★★★★★♦★★★★★★★It************************** 


common/ numrec/  nf us , nf ss , ndsw, ncs t , nss , ness , nudg , nvdg , nu l w , nrs 
*,  ncl t ,ncls,ntsf ,nqf l , ncum.nl sp,nps,nrma,nt gf ,nrmg,ntcc,nth2o 
*,  nevp,nctt,ncts,ncsa,ncla,ncta,npcp,nvts,nvtg,npme 
*,  nsno,nwet,nwev,nalb,nblk(4) 

*,  ntdiab.nqdi ab,nudiab,nvdi ab,ntcum,nqcum,nucum,nvcum 
* ,  ncumf , nugwd , nvgwd, nt  radh , nt l wr , nt  radc , nt l wre 
*,  navcld,navcvc,ndthd,ndqhd,nduhd,ndvhd 

*,  nslp,nslpv,ntmp,ntmpv,nhgt,nhgtv,nsht,nshtv,nuwn,nuwnv 
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*,  nvun,nvunv,npsi , npsiv.nchi ,nchiv 

*,  ntmpcs,nshtcs,nrhcs,ncldcs,nuucs,nvvcs,nstmcs 

dimension  nrecs(numave) 
equivalence  (nrecs.nfus) 

common/ rec lab/  Ifws, Ifss, Idsw, lest, Iss, less, ludg, Ivdg, lulu, Irs 
*,  Iclt, Icls, Itsf ,  Iqf  l, l cum,  lisp,  tps,  Irma,  Itgf ,  Irmg,  Itcc,  lth2o 
*,  levp, lett, lets, lesa, Icla, Icta, Ipcp, Ivts, Ivtg, Ipme 
*,  Isno, luet, luev, latb, lbtk(4) 

*,  l tdiab( (pout), lqdiab( (pout), ludiab< tpout) , l vdiab( l pout) 

* ,  l tcum( l pout ) , lqcum( l pout ) , l ucum{ l pout ) , l vcum( l pout ) 

*,  lcumf(lpout),lugwd(lpout),  lvgud(lpout),ltswh(lpout) 

*,  ltlur( l pout), l tswc< l pout), Itlwrc(lpout) 

*,  lavcld(lpout), lavcvc( Ipout), IdthdC Ipout), ldqhd( Ipout ) 

*,  lduhd( Ipout), ldvhd( Ipout ) 

*,  l sip, Islpv, l tmpC Ipout), ItmpvC Ipout), lhgt< Ipout), Ihgtv(lpout) 
*,  lsht( Ipout), Ishtv(lpout), luun( Ipout ), luwnvC Ipout) 

*,  lvwn( Ipout), l vwnv( Ipout), IpsiC Ipout ), IpsivC Ipout ) 

*,  lchi( Ipout), IchivC Ipout) 

*,  Itmpcs,  Ishtcs,  Irhcs,  l cldcs,  luucs,  Iwcs,  Istmcs 

character*8  Ifus, Ifss, Idsw, lest, Iss, less, ludg, Ivdg, lulu, Irs 
*,  Iclt,  Icls,  Itsf,  Iqf  l ,  leum,  l  Isp,  Ips,  Irma,  Itgf ,  l ring,  Itcc,  lth2o 
*,  levp, lett, lets, lesa, Icla, Icta, Ipcp, Ivts, Ivtg, Ipme, tsno, luet 
*,  luev, lalb, Iblk 

*,  ltdiab, Iqdiab, ludiab, Ivdiab, Itcum, Iqcum, lucum, Iveum 
*,  leumf , lugwd, Ivgud, Itsuh, Itlur, Itsuc, Itlurc 
*,  lavcld, lave vc, Idthd, Idqhd, Iduhd, Idvhd 

*,  Islp, Islpv, Itmp, Itmpv, Ihgt, Ihgtv, Isht, Ishtv, luwn, luunv 
*,  Ivwn, Ivunv, Ipsi , Ipsi v, Ichi , Ichiv 

*,  Itmpcs,  Ishtcs,  Irhcs,  Icldcs,  luucs,  Iwcs,  Istmcs 

character*8  drecs(numave),clred(21*lpout),dreca(2+K*lpout) 
*,  clcros(7) 

equivalence  (clrecs, Ifus) 
equivalence  (cl reel , ltdiab) 
equivalence  (clreca,lslp) 
equivalence  ( cl cros, Itmpcs) 

dimension  pout ( Ipout ),psout( Ipout) , pcs l ( Im) 

logical  superg.crf , fpass 

character*64  capsun 
character*36  maplst 
character*8  i dent (24) 
character*8  idtgs, idtge, idtgx, idx 
character*8  l rec 
character*^  iplev 

character*48  hccn,c3ave,f i larg.fi Isml , trpf i l ,c3grid, shist 

character*48  c3path,avpath 

character*24  llab(numave) 

character*2  month, neumon 

character*6  c3av,shst 

character‘7  hcc,c3g,trpf 

character*64  cap(7) 

dimension  plotc(5,numave) 
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data  zero  /0.0/ 

data  pout/10.0,100.0,200.0,500.0,850.0.1000.0/ 
data  pcs/1 0.0, 20. 0,30. 0,50. 0,70.0, 100.0, 150.0,200.0, 250.0 
*,  300.0,400.0,500.0,600.0,700.0,800.0,900.0,950.0,1000.0/ 
c 

open(2,fi  ^e=,/u/b/rosmond/timave/labpa^^l,  ,form='  formatted' ) 
c 

read(2,800)  (clrecs(k), l lab(k), (plotc(ki ,k),ki=1,5) 

*,  k=1,numave) 

800  format(6x,a8,1x,a24,2x,4f8.3,f4.1) 
c 

print  850,(clrecs(k), llab(k) 

*,(plotc(ki ,k),ki=1,5),k=1 .nunave) 

850  format(6x,a8,1x,a24,2x,4f14.5,f4.1) 
c 

isize=  -1 

c 

do  2  k=1,  nunave 
nrecs(k)=  k 

2  continue 
c 

do  3  k=1 , Ipout 
psout(k)=  log(pout(k)) 

3  continue 
c 

do  4  k=1 , Im 
pcsl(k)=  log(pcs(k)) 

4  continue 
c 

c  compute  interpolation  coeffs  and  tranpose  array  y 
c 

pi=  4.0*atan(1 .0) 
pio2=  pi *0.5 
rad=  6375000. 
radsq*  rad*rad 
f  pass*. true, 
c 

call  fftfax(im,ifax, trigs) 
c 

call  gausl3( jm,-1 .0,1 .0,weight,sinl ) 

do  5  j=2,jm*1 

alat(j)=  asin(sinl( j-1 )) 

5  continue 

tem=  pio2*1.0e-6 
alat(1)=  -tern 
alat(jm+2)=  tern 

cal l  sortml( jtrun,mlmax,msort, l sort, ml  sort) 
c 

do  152  ml=1,mlmax 
rl=  Isort(ml) 
rm=  msort(ml)-1 
rlm=  rl-1.0 

ifOnsort(ml).eq.l)  rm=  0.0 
ifdsort(ml).eq.l)  rlm=  0.0 
eps4(ml)=  rl*rlm/radsq 
cim(ml)=  rm 
152  continue 
COIRS  IVDEP 

do  155  j=1,jm/2 

sinl( jm+1 - j)=  -sinl(j) 

ocos(j)=  1.0/( 1. 0-si nl(j)*s ini ( j)) 

ocos( jm*1 - j )=  ocos(j) 

cosl(j)=  1.0/sqrt(ocos(j)) 

cosl( jm*1- j)=  cosl(j) 

do  155  i=1,im 

onocos(i,j)=  ocos(j) 
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onocosd,  jm-j*1)=  ocos(j) 
155  continue 


call  lgndr( jm/2, jtrun, idim2,mlsort,poly,dpoly,sinl ) 


do  12  k=1,24 
ident(k)=‘  1 

12  continue 

opend.f  ile='pcards',fonn='  format  ted1 ) 

read<1, 80)  filarg 
readd, 80)  distill 
read(1, 80)  avpath 
readd, 80)  c3path 
80  format(a48) 

call  chlen(avpath, Ipath) 

call  chlen(c3path,l3path) 

c3av='/c3ave' 

hcc=  '/hccn32‘ 

trpf='/trpf i l ' 

c3g=  '/c3grid' 

shst='/shist' 

call  cfopen(fi larg, 41760, istat) 
call  cfopenCfi Isml, 10512, istat) 

read(1, 100)  capsun 
100  format (a64) 

print*,'  capsun  ', capsun 

120  readd,  1 40, end=  160)  maplst 
140  format(a36) 

print*,'  maplst  '.maplst 
iszold=  isize 

if  (maplstd  :8).eq. 'nomodata' )  go  to  160 
read(maplst,150)  tree, isize, idtgs, idtge 
150  format (a8, lx, i 1 , 1x,a8, 1x,a8) 
superg=  isize. gt.O 
idtg(  1  )= '  7901 01 00 ' 

find  itaus  -..id  i^ue  corresponding  to  beginning  and  ending 
dtg  of  aver  ging  period 

print*,  '  idtg  '.idtgd) 
print*,  '  idtgs  '.idtgs 
print*,  '  idtge '.idtge 
call  tothrs(idtg(1),itaux) 
call  tothrsddtgs,  itaus) 
call  tothrsl idtge, itaue) 
print*,'  taux  '.itaux 
print*,'  taus  '.itaus 
print*,'  taue  '.itaue 
itaus=  itaus-itaux+24 
i taue=  i taue- i taux+24 
print*,'  taus  '.itaus 
print*,'  taue  '.itaue 

month=idtgs(3:4) 

c3ave=  c3path(1 : l3path)//month//c3av 
hccn=  avpathd :  lpath)//month//hcc 
trpf  i  1=  avpathd :  lpath)//month//trpf 
c3grid=  avpathd :  lpath)//month//c3g 
shist*  avpathd : lpath)//month//shst 


c3ave=  c3path(1 : I3path)//c3av 
hccrv=  avpathd :  lpath)//hcc 
trpfil*  avpath(1 : lpath)//trpf 
c3grid=  avpathd :  Ipath)/ 
shist*  avpathd : l path )//shst 
c 

idtgx=  idtgs 
c 

i tx=  -1 

call  nfreadChccn.1,2,1,1, Iccn.c.itx, idx, is tat) 
c 

if(fpass)  then 
fpasss. false, 
c 

call  nfread(trpfil,1,2,1,3, imjm.pdiff ,-1,idx, istat) 
call  nfread(c3grid,1,2,1,1,imjm,sgeo,itaus,idx,istat) 
call  qpnhCpdiff , 'pdiffpdiff ‘,1,1,144,72) 
c 

call  zi 1 chCphi ps, imjm*( lpout*7+1 )*2) 
call  zi lch(ave2d, imjm*36) 
call  zi lch(ave3d, imjm*lpout*21 ) 
call  zi lch(tempcs, imjm*lm*7) 
do  351  i=1,imjm 
dswi i i C i )=  0.0001 
olriii(i)=  0.0001 
351  continue 
c 

ipt=  0 
numav=  0 

do  165  itau=  itaus, itaue,24 
itau1=  itau-18 
itau2=  itau-12 
itau3=  itau-6 
itau4=  itau 
c 

c  spectral  fields 
c 

call  nf readt sh i st , i pt , 2 , 1 , nspec , i di m2 , x 1 1 , i t aul , i dx , i st at ) 
call  nfread(shist, ipt,2, 1 .nspec, idim2,xl2, itau2, idx, istat ) 
call  nfreadCshist, ipt,2, 1 , nspec, idim2,xl3, i tau3, idx, istat) 
cal l  nfread(shist, ipt,2, 1 .nspec, idim2,xl4, itau4, idx,  istat) 
c 

do  110  i=1, idim2*nspec 

vornowC i ,1 , 1 )=  0.25*(xl1(i , 1 )+xl2< i ,1)+xl3( i , 1 )+xt4< i , 1 )) 
110  continue 
c 

call  transr(lm,poly,vornoH,rvor) 
call  transr(lm,poly,divnoH,rdiv) 
call  tranuvC Im, radsq,ocos,eps4,cim,poly,dpoly 
*,  vornoH,divnoH,ut,vt) 
c 

do  115  k=1 , Im 
vornow(1,1,k)=  0.0 
vornoMd  ,2,k)=  0.0 
divnowd,1,k)=  0.0 
divnoH(1,2,k)=  0.0 
do  115  ml=2,mlmax 

vornowCml , 1 , k)»* 1 .0e-6*vornow(ml , 1 , k)/eps4<ml ) 
divnoM(ml,1,k)=-1.0e-6*divnow<ml,1,k)/eps4(ml) 
vornoM(ml,2,k)=-1 .0e-6*vornoH(ml,2,k)/eps4(ml ) 
divnowCml ,2,k)=-1 .0e-6*divnow(ml ,2,k)/eps4(ml ) 

115  continue 

call  transr(lm,poly,temnow,tt) 
call  transr(lm,poly,shnow,sht) 
call  transrd , poly, plnoH.pt) 
call  transr(lm,poly,vornow,rstrm) 
call  transr(lm,poly,divnow,vpot) 
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cat  t  trngra(1,cim, poly, dpoly.pt non, dipt ,dtpt ) 
c 

call  p2kapatimjm, lm,asig,bsig,pt,pk,plog,pk2,pt t ) 
c 

c  hydrostatic  equation 
c 

cp=  1005.45 
pp=  log(IOOO.O) 
do  402  i =1 , imjm 

phi  t i , Im)  =  sgeol i )+cp*tt( i , lm)*tpk2t i , lm)-pkt i , Im) ) 
sht( i , lm)=  sht( i , lm)/tdasi gt tm)+dbs i gt lm)*pt t i , 1 ) ) 

402  continue 

do  405  1=1, la- 1 
do  405  i=1 , imjm 

sht(i,l)=  sht( i , l )/(dasig( l )+dbsigf l )*ptf i ,  1 )) 
phi(i,l)=  tt(i,l+1)*(pk(i, 1+1 )-pk2( i , l))+tt( i , l ) 

*  *(pk2( i , l )-pk( i ,  l ) ) 

405  continue 

do  440  l rev  =  1 , lm-1 
l  =  Im  -  l rev 
do  440  i=1 , imjm 

440  phi ( i , l )=phi ( i , 1+1 )+phi ( i , l )*cp 
c 

cal l  womega(im, jm, Im, onocos.bsig, das ig, dbsig.pt 
*,  dtpl,dlpl,rdiv,ut,vt,xl) 
c 

c  compute  specific  humidity  and  temperature 
c 

do  330  k=1,lm 

do  326  i=1 , imjm 

sht(i,k)=  maxtshtti,k),1.0e-6) 

tt(i ,k)=  tt< i ,k)*pk( i ,k)/<1 .0*0.608* sht(i ,k)) 

326  continue 
c 

call  qsatvs( imjm,tt(1,k),pll(1,k),x) 
c 

do  330  i=1 , imjm 

rhcs( i ,k)=  rhcst i ,k)+mint 1 .O.shtt i ,k)/xt i >) 
tempcs(i,k)=  tempest i ,k)+tt(i ,k) 
shtcs(i,k)=  shtcsC i ,k)+sht(i ,k) 
uucs(i,k)=  uucs(i,k)+ut(i,k) 
wcs(i,k)=  west  i  ,k)+vtt  i  ,k) 
strmcs(i,k)=  strmcsl i ,k)+xt  ( i  ,k) 
sht(i,k)=  logtshtci , k) ) 

330  continue 
c 

kmid=  4 

call  svarfimjm, lm,kmid,pt,sgeo,pdiff , t save, phi , tt.tmppl ,slpx) 
call  geotrpl im, jm, Im, Ipout.psout, tt.phi ,sgeo,slpx 
*,  plog.xl.xll) 
c 

call  avevartimjm, Ipout.xl, temps, tempsv) 
call  avevartimjm, Ipout.xl 1, phi ps, phi psv) 
call  avevartimjm, l.slpx, sip, slpv) 
c 

call  sZptrpt im, jm, Im, tpout ,ut,utt 1 , lm),plog,xl ,psout,0.5) 
ca 1 1  avevar  t i m j  m, l pout , x l , upres , upres v) 
cal l  s2ptrpf im, jm, Im, Ipout.vt, vtf 1, lm),ptog,xl,psout,0.5) 
call  avevartimjm, lpout,xl,vpres,vpresv) 
cal l  s2ptrpt im, jm, tm, lpout,sht,sht<1, tm),plog,xl ,psout,0.0) 
do  333  k=1,lpout 
do  333  i=1 , imjm 
xl(i,k)=  exptxtti,k))*1000.0 
333  continue 

cal l  avevar t imjm, Ipout.xl ,sphum,sphumv) 

call  s2ptrptim, jm,lm, l pout , rstrm, rstrmt 1 , Im) .plog.xl ,psout ,0.5 ) 
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call  avevarfimjm, lpout,xl,psi,psiv) 

call  s2ptrp(im, jm, Im, lpout,vpot,vpot(1, lm),plog,xl,psout,0.5) 
call  avevarfimjm, lpout,xl,chi ,chiv) 
c 

c  read  ground  wetness  and  snow 
c 

call  nf read(c3grid, ipt,2,6,1, imjm, x 1(1 , 1 ), itaul , idx, istat) 
call  nfread(c3grid,ipt,2,6,1, imjm,xl(1 ,2), itau2, idx, istat) 
cal l  nfread(c3grid, ipt,2,6, 1, imjm,xl(1,3),i tau3, idx, istat) 
call  nfread(c3grid, ipt,2,6,1, imjm,xl(1 ,4), i tau4, idx, istat) 
c 

call  nf read(c3gr id,  ipt, 2,9,1 ,  imjm, x 1(1  ,5),  itaul  ,idx,  istat) 
call  nfread(c3grid, ipt, 2,9,1, imjm, xl(1 ,6), itau2, idx, istat) 
call  nfread(c3grid, ipt, 2, 9, 1, imjm,xl(1 ,7), itau3, idx, istat) 
call  nfread(c3grid, ipt, 2, 9, 1, imjm.xl ( 1 ,8), i tau4, idx, istat) 
c 

do  337  i=1,imjm 

xl(i,nwet)=0.25*(xl(i,1)+xl(i ,2)+xl(i,3)+xl(i,4)) 
ave2d(  i ,  nwet  )=ave2d(  i ,  nwet  )+x  l  { i ,  nwet ) 

ave2d(i,nsno)=ave2d(i,nsno)+0.25*(xl(i,5)+xl(i,6)+xl(i,7)+xl(i,8)) 
337  continue 
c 

c  read  2-d  daily  mean  history 
c 

cal l  nfread(c3ave,1 ,2, 1,20, imjm.xl , i tau, idx, istat) 
c 

cal l  solar(xl(1 ,nalb), idtgx,sinl,cosl) 
c 

do  340  i=1,imjm 

ave2d( i ,nf ws)=  ave2d( i ,nf ws)+xl ( i ,nf ws) 
ave2d( i ,nfss)=  ave2d( i ,nfss)+xl(i ,nfss) 
ave2d( i , ndsw)=  ave2d( i , ndsw)+x l < i , ndsw) 
a ve2d( i , nss ) =  a ve2d( i , nss )+x l ( i , nss ) 
a ve2d( i , nudg ) =  a ve2d( i , nudg ) +x l ( i , nudg ) 
ave2d( i ,nvdg)=  ave2d( i ,nvdg)+xl ( i ,nvdg) 
ave2d( i ,nulw)=  ave2d( i ,nulw)+xl( i ,nulw) 
ave2d(i ,nrs)=  ave2d( i ,nrs)+xl( i ,nrs) 
ave2d( i , ntsf )=  ave2d( i , ntsf )+x l ( i , ntsf ) 
a ve2d( i , nqf l )=  ave2d( i , nqf l )+x l ( i , nqf l ) 
ave2d( i ,ncun)=  ave2d( i ,ncum)+xl(i ,ncum) 
ave2d( i ,nlsp)=  ave2d( i ,nlsp)+xl ( i ,nlsp) 
a ve2d(  i ,  nps  )=  ave2d(  i ,  nps  )*x  l  ( i ,  nps ) 
ave2d( i , nrma )=  ave2d< i , nrma )+x l ( i , nrma ) 
ave2d( i , ntgf )=  ave2d( i , ntgf )+x l ( i , ntgf ) 
ave2d( i ,nrmg)=  ave2d( i ,nrmg)+xl( i ,nrmg) 
ave2d( i ,nalb)=  ave2d( i ,nalb)+xl(i ,nalb) 

340  continue 
c 

c  compute  mean  square  ground  and  surface  air  temps  for  later 
c  variance  determination;  also  ground  wetness 
c 

do  345  i=1 , imjm 

ave2d( i ,nvts)=  ave2d( i ,nvts)+(xl( i ,ntsf )-273. 15)**2 
ave2d(i ,nvtg)=  ave2d(i ,nvtg)+(xl(i ,ntgf)-273.15)**2 
ave2d( i ,nwev)=  ave2d( i ,nwev)+xl ( i ,nwet)**2 
345  continue 
c 

c  read  3-d  daily  mean  history 
c 

do  350  k=1 ,21 
kk=  21+<k-1)*lm 

call  nf read(c3ave, ipt,2,kk, Im, imjm,tt,itau, idx, istat) 
c 

if(k.eq.9)  then 
do  352  i=1 , imjm 

ave2d( i ,ntcc)=  ave2d( i ,ntcc)+tt( i , 1 ) 
ave2d( i ,nth2o)=ave2d( i ,nth2o)+tt( i ,2) 
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352  continue 

call  zilch(tt,imjm*2) 
endif 
c 

if(k.eq.7)  then 
c 

c  compute  cloud  forcing  using  method  3 
c 

do  353  i=1,imjm 
erfiii(i)*  tt(i,1) 
if(crfiii(i).gt.1.9)  then 

ave2d(i ,ncst)=  ave2d( i ,ncst)+xl( i ,ndsw)-xld ,ncst) 
ave2dd  ,ncss)=  ave2d( i ,ncss)+xl(i ,nss)-xl( i ,ncss) 
dswi i i ( i )=  dswi i i ( i )♦! .0 
endif 
c 

c  l ongwave 
c 

if (erf i ii(i ).gt. 0.999)  then 

ave2d( i , nc 1 1 ) =  ave2d( i , nc 1 1 )+x  U i , nc 1 1 ) - x l ( i , nul  w ) 
ave2d(i  ,ncls)=  ave2d(i  ,ncls)+xl(i  .nclsJ-xK  i  ,nrs) 
olriii(i)=  olriii(i)+1.0 
endif 

tt(i,1)»  0.0 

353  continue 
endif 

c 

if(k.eq.16)  then 
do  357  1=  1, Im 
do  357  i=1,imjm 
xl(i,l)=  tt(i,l) 

357  continue 
c 

else  if'*.eq.17)  then 
c 

do  359  1=1,  Im 
do  359  i=1 , imjm 

cldscs(i ,  l  )  =  cldscs< i , l )+max(xl( i , l),tt(i , l)) 

359  continue 
c 

c  else  if(k.eq.15)  then 
endi  f 
c 

call  f ixdat(k,tt,tt, imjm, Im) 
c 
c 

call  s2ptrp( im, jm, Im, l pout , tt , tt(1 , Im) ,plog, rvor,psout ,0.5) 
c 

do  355  j=1 , l pout 
jj=  j+(k-1 )*lpout 
do  355  i=1 , imjm 

ave3d( i , j  j )=  ave3d( i , j  j )+rvor( i , j ) 

355  continue 
c 

350  continue 
c 

numav=  numav+1 
c 

call  dtgf ixlidtgx, idtgx,24) 
newmon=idtgx(3:4) 
c 

if (newmon.ne. month)  then 
month=  newmon 

c3ave=  c3path(1 : l3path)//month//c3av 
hccn=  avpath( 1 : Ipath )//month//hcc 
trpf  i  1=  avpathd:  lpath)//month//trpf 
c3grid=  avpathd :  Ipath  )//month//c3g 
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shist=  avpath(1 : lpath)//month//shst 
endif 

165  continue 

fac=  1.0/numav 
do  167  i =1 ,  imjrn 

ave2d( i ,nfws)=  ave2d(i ,nfus)*fac 
ave2d( i ,nfss)=  ave2d( i ,nfss)*fac 
ave2d( i ,ndsw)*  ave2d( i ,ndsw)*fac 
ave2d( i ,ncst)=  ave2d( i , nest )/dswi i i ( i ) 
ave2d( i ,nss)=  ave2d(i ,nss)*fac 
ave2d( i ,ncss)=  ave2d( i ,ncss)/dswi i i ( i ) 
ave2d( i , nudg)=  ave2d( i ,nudg)*fac 
ave2d( i ,nvdg)=  ave2d( i ,nvdg)*fac 
ave2d(i,nulw)=  ave2d(i,nutw)*fac 
ave2d( i ,nrs)=  ave2d(i,nrs)*fac 
ave2d( i ,nclt)=  ave2d( i ,nclt)/olri i i ( i ) 
ave2d( i ,ncls)=  ave2d( i ,ncls)/olri i i ( i ) 
ave2d( i ,ntsf )=  ave2d( i ,ntsf )*fac 
ave2d( i ,nqf l )=  ave2d( i ,nqf t )*fac 
ave2d(i,ncum)=  ave2d(i,ncum)*fac 
ave2d( i ,nlsp)=  ave2d( i ,nlsp)*f ac 
ave2d( i ,nps)=  ave2d( i ,nps)*fac 
ave2d(i,nrma)=  ave2d(i,nrma)*fac 
ave2d( i , ntgf )=  ave2d( i , ntgf )*f ac 
ave2d( i ,nrmg)=  ave2d(i,nrmg)*fac 
dswi i i ( i )=  dswi i i ( i )*f ac 
olriii(i)=  olriii(i)*fac 

167  continue 

temx=  3600.0*24. 0/2. 52e6 
do  260  i =1 , imjm 

ave2d( i , nett >=  ave2d( i , nc l t)+ave2d( i , nest) 

ave2d( i ,ncts)=  ave2d( i ,ncls)+ave2d( i ,ncss) 

ave2d( i , nesa ) =  ave2d( i , nest ) - ave2d( i , ness ) 

ave2d(i ,ncla)=  ave2d(i ,nclt)-ave2d(i ,ncls) 

ave2d( i , ncta) =  a ve2d( i , nesa )+ave2d( i , nc l a) 

ave2d(  i  ,npcp)=  ave2d(  i  ,nciin)+ave2d(  i  ,nlsp) 

ave2d(i,nevp)=  temx*ave2d(i ,nfws) 

ave2d( i , ntsf )=  ave2d( i , ntsf ) - 273 . 1 5 

ave2d( i , ntgf )=  ave2d( i , ntgf ) - 273 . 1 5 

ave2d( i , npme)=  a ve2d( i , npep) - ave2d( i , nevp) 

ave2d(i ,nvts)=  fac*ave2d( i ,nvts)-ave2d( i ,ntsf )**2 

ave2d( i ,nvtg)=  fac*ave2d( i ,nvtg)-ave2d(i ,ntgf )**2 

ave2d( i ,nsno)=  ave2d( i ,nsno)*fac 

ave2d( i , nuet )=  ave2d( i ,nwet ) *  f ac 

ave2d( i ,nwev)=  fac*ave2d(i ,nwev)-ave2d(i ,nwet)**2 

ave2d( i ,ntcc)=  min(0.999,ave2d( i ,ntcc)*fac) 

ave2d( i , nth2o)=ave2d( i ,nth2o)*f ac 

ave2d{ i ,nalb)=  ave2d(i ,nalb)*fac 

if (ave2d( i ,nalb).gt .0.0) 

*  ave2d( i ,natb)=  1.0-ave2d(i,ndsw)/ave2d(i,nalb) 
tmpp1(i)=  tmppl ( i )*nunav 

260  continue 

endif 

call  s2pcs( im, jm, lm, tenpes, tmppl ,plog,cstt ,pcsl , fac) 

call  s2pcs(im, jm, lm,shtcs,shtcs(1, lm),plog,cssht,pcsl,fac) 

call  s2pcs(im, jm, lm,rhcs,rhcs(1, lm),plog,csrh,pcsl,fac) 

cal l  s2pcs( im, jm, lm,cldscs,cldscs(1, lm), pi og, esc  Ids, pcs l ,fac) 

cal l  s2pcs( im, jm, lm,uucs,uucs( 1 , lm),plog,csuu,pcsl , fac) 

call  s2pcs(im,  jm,  lm,wcs,wcs(1,  lm),plog,csvv,pcsl,fac) 

call  s2pcs(im, jm, lm,strmcs,strmcs(1, lm),plog,csstrm,pcsl ,fac) 
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facx=  1.0/im 
do  262  j=1,jm 
ptzm(j)=  0.0 
do  264  i=1,im 
ptzm(j)=  ptzm(j)+pt(i,j) 

264  continue 

ptzm(j)=  ptzm(j)*facx 
262  continue 
c 

cal l  hadly(im, jm, Im, jtrun,mlmax, weight, poly, msort.csstrm) 
c 

do  265  k=1 , Im 
do  265  j=1,jm 

cscldsC j,k)  =  max(0.0,csclds( j,k)) 
csrh(j,k)=  max(0.0,min(1 .0,csrh( j,k))) 
csuu( j,k)=  csuu(j,k)*rad/cosl(j> 
esw(j,k)=  csw(j,k)*rad/cosl(j) 

265  continue 
c 

if (Irec.eq. 'CROSSECS' )  then 
do  362  i=1,7 

cap( i )(1 :24)=  llab( i+numave-7) 
cap(i) (25:64)=' 

362  continue 

capsun='AMJP  TIMEAVE  FROM  '//idtgs//'  TO  '//idtge//'  ' 

lencs=  7*lm*jm+lm+ jm 
read( idtgs, 1 (2i2) 1  )i1 ,  i2 
itau=  100*i2+i1 

call  dgwri t( jm,7, lencs, i tau, idtgs, capsun, cap, cstt, 'zmeancss ‘ > 
go  to  120 
end  if 
c 

c  ****************************************************************** 

c 

do  15  kid=1,36 
if(lrec.eq.clrecs(kid))  then 
do  17  i=1,imjm 
avsum(i)=  ave2d(i,kid) 

17  continue 
c 

kl=  kid 
klev=1 

iplev='  1 
go  to  52 
end  if 

15  continue 
c 

do  57  kid=1 ,21*lpout 
ifOrec.eq.clrecl(kid))  then 
kl=  40+kid 
iplev=  lrec(5:8) 
read(iplev, '(i4)')  ip 

if(ip.gt.O)  call  levget(ip,pout,lpcut,klev) 
print*,  kid, Irec.klev 
do  59  i=1,imjm 
avsum(i)=  ave3d(i ,kid)*fac 
59  continue 
go  to  52 
endif 

57  continue 
c 

do  55  kid=1,2+14*lpout 
c 

if(lrec.eq.clrecaCkid))  then 
c 

kl*  40+21*lpout+kid 
iplev=  lrec(5:8) 
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readd  pi  ev,  *  <  i 4 ) ' )  ip 

if(ip.gt.O)  call  levget(ip,pout,  Ipout.k  ...-v) 
print*,  kid, Irec.klev 

sea  level  pressure 

if(lrec.eq.lslp)  then 
call  ave(imjm,fac, sip, avsum) 
klev=  1 


else  if ( Irec.eq. Islpv)  then 

cal l  var< imjm,fac,slp,slpv, avsum) 

klev=  1 

temperature 

else  ifdrec.eq. Itmp(klev))  then 

ca 1 1  ave( i m j m, f ac , temps( 1 , k l ev) , avsum) 

else  ifdrec.eq.  Itmpv(klev))  then 

call  var(imjm,fac,temps(1,klev),tempsv(1,k lev), avsum) 

geopotential 

else  if ( Irec.eq. Ihgt(klev))  then 
call  ave(imjm,fac,phips(1,k lev), avsum) 

else  ifdrec.eq.  Ihgtv(klev))  then 

call  var ( im jm, f ac , ph i ps( 1 , k l ev) , ph i ps v( 1 , k l ev ) , avsum) 

specific  humidity 

else  if(lrec.eq.lsht(klev))  then 
call  ave(imjm,fac,sphum(1,k lev), avsum) 

else  ifdrec.eq.  Ishtv(klev))  then 

call  var( imjm,fac,sphum(1,k lev), sphumv(1,k lev), avsum) 

u  component  Kind 

else  ifdrec.eq.  luwn(klev))  then 
call  ave( i m j m, f ac , upresC 1 , k l ev) , avsum) 


else  ifdrec.eq. luwnv(klev))  then 

call  var(imjm,fac,upres(1,klev),upresv(1 ,k lev), avsum) 
do  210  i=1 , imjm 

avsun(i)=  avsumf i )*radsq*ocos( i ) 

210  continue 

v  component  wind 


else  ifdrec.eq.  Ivwn(klev))  then 

cal l  ave( imjm, f ac , vpres( 1 , k lev) , avsum) 


else  ifdrec.eq.  Ivwnv(klev))  then 

cal l  var( imjm, f ac, vpres( 1 , k l ev) , vpresv( 1 , k l ev) , avsum) 
do  220  i=1,imjm 

avsum( i )=  avsum( i )*radsq*ocos( i ) 

220  continue 

stream  function 


else  ifdrec.eq. Ipsi(klev))  then 
cal l  ave(imjm,fac,psi(1 ,k lev), avsum) 

else  if (Irec.eq. Ipsiv(klev))  then 

call  var ( im j m, f ac , ps i ( 1 , k l ev ) , ps i v( 1 , k l ev ) , avsum) 
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c  velocity  potential 
c 

else  if(lrec.eq.tchi(klev))  then 
call  ave( imjm, f ac , ch i ( 1 , k l ev) , avsum) 
c 

else  if ( Irec.eq. lchiv(klev))  then 
cal l  var< imjm, f ac,chi ( 1, k lev), chi v(1 ,klev), avsum) 
c 

endif 
go  to  52 
endi  f 

55  continue 
c 

go  to  120 
c 

52  continue 
c 

call  sumcosC im, jm, cos l , avsum ,gmean) 
print*,'  gmean=  ' ,gmean 
c 

ci=  plotc(1,kl) 
co=  plotc(2,kl) 
add=plotc(3,kl ) 
amp=plotc(4,kl ) 
c 

call  lablsCkl, lrec(1 :3), iplev, isize, idtg(1 ), idtgs, idtge.gmean 
* ,  capsunl 1 : 24 ) , c i , co, add, amp, l l ab( k l ) , i dent ) 
print*,  ident 
c 

if(superg)  then 
imax=  288 
jmax=  145 
else 

imax=  144 
jmax=  73 
endif 
c 

ik=  1 

if(isize.ne.iszold)  ik=  0 

call  hterpC ik, avsum, im, jm, glob, x in, yin, imax, jmax, a l at , vwk) 
if(lrec(1:3).eq.luwn(klev)(1:3))  call  wndfixlglob, imax, jmax) 
if(lrec(1:3).eq.lvwn(klev)(1:3))  call  wndfixCglob, imax, jmax) 
if(plotc(5,kl).eq.1)  call  f ixi tlglob, imax* jmax, 0.0) 
cal l  qprnthfglob, ident, 1 , 1 , imax, jmax) 
call  qpnh(glob, ident, 1 , 1 , imax, jmax) 
c 

lenx=  imax* jmax 
if(superg)  then 

c  call  ritasc(ident, lenx.glob) 

cal l  cfwri t(f i larg, ident, glob, lenx,0, i stat ) 
else 

call  gloflpCimax, jmax, glob, gf lip) 
c  call  ri tasc( ident, lenx,gf l ip) 

cal l  cfwri t(f i Isml , ident,gf l ip, lenx,0, istat ) 
endi  f 
c 

go  to  120 
c 

160  continue 
stop 
end 

subroutine  f ixdat(kk,tt,mtt, imjm, Im) 
c 

dimension  mttfimjm, lm),tt(imjm, Im) 
c 

i count =  0 
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c 


do  10  k=1 , Im 
do  10  i=1,in)jm 
ll=ibits(mtt(i,k),0,24) 
if(ll.ne.O)  then 
tt(i,k)=  0.0 
i count =  icount+1 
end  if 

10  continue 

if(icount.eq.O)  return 

print*, kk,'  fixdat  found  '.icount,'  bad  values' 

return 

end 

subroutine  ave(imjm,fac,x,avsum) 
c 

dimension  x(imjm),avsi«n(imjm) 
c 

do  10  i=1 , imjm 
avsun(i)=  fac*x(i) 

10  continue 
return 
end 

subrout i ne  avevar ( i m j  m, 1 1 , x l , ave , va  r ) 
c 

dimension  ave(imjm,  1 1 ) ,  var(  imjm,  ll),xl(imjm,  U) 
c 

do  10  k=1,ll 
do  10  i=1,imjm 
ave( i  ,k)=  ave( i ,k)+xl( i ,k) 
var(i,k)=  var(i,k)+xl(i ,k)**2 
10  continue 
c 

return 

end 

subroutine  bcubv< iderv.f , ix, jy,dout,mn,yr,yin, ten) 
c 

c  parameter  list 
c 

c  f , ix, jy,dout,mn,xin,yin  -  see  below 
c 

dimension  f (ix, jy),dout(mn),yr( ix, jy).yin(mn) 
c 

dimension  fxx( ix, jy),fyy( ix, jy)  pjy(mn,4), tpl (mn,4) 
dimension  ipt(mn), jpt(mn) 
c 
c 

c  a  bicubic  spline  interpolator  to  interpolate  from  a  grid 

c  with  constant  i  (first  dimension)  grid  spacing  and  variable 

c  j  (second  dimension)  grid  spacing  tO  a  grid  with 

c  variable  grid  spacing,  all  grids  are  assumed  to  have  point 

c  (i,1)  in  the  lower  left  corner  with  i  increasing  to  the  right 

c  and  j  increasing  upward, 

c 
c 

c  compute  ipt.jpt  and  pjy 

c 

call  setupv( iderv,yr,yin,mn, ix, jy.pjy, ipt.jpt, tpl) 

jm=  mn/ix 

ixm1=ix-1 

jym1=jy-1 

jym2=jy-2 

i jm2=ix*jy-2 

ixjym2=ix*jym2 

mn4=mn*4 


c 

c 

c 


compute  fyy 
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11=  ixjym2+ix 
do  100  i=1 , 1 1 
fxx(i,2)=  yr(i,2)-yr(i,1) 

100  continue 
c 

do  210  i=1,ixjym2 

fyy(i,2)=  ten*<fxx(i,3)*(f(i,1)-f(i,2))+fxx(i,2)*(f(i,3)-f(i,2))) 
*  /(fxx(i,3)*fxx(i,3)*fxx(i,2)) 
fxx(i,2)=  fxx(i,2)/fxx(i,3) 

210  continue 
c 

call  trdiwlix, jym2,fxx(1 ,2),fyy(1 ,2)) 
c 

call  zi lch(fyy, ix) 
call  zi lch(fyy(1 , jy) , ix) 
c 

call  gathv(mn, ix, jy, ipt, jpt ,fyy,f ,tp1 ) 
c 

do  130  i=1,mr>4 

tpl ( i , 1 )=tp1 ( i , 1 )*p jy( i , 1 ) 

130  continue 

do  140  i=1,mn 

dout(i)=tp1<i,1)+tp1(i,2)+tp1(i,3)+tp1(i,4) 

140  continue 
return 
end 

subroutine  dgwri t( idl , id2, len, i tau,cdtg,capsun,cap,x, label ) 
c 

dimension  x( len) 
character*8  cdtg 
character*8  label 
character*48  filedg 
character*14  cform 
character *64  cap(id2),capsun 
C 

character*32  dagdir 
c 

dagdi r='/u/b/rosmond/pcmdi/pcmf Ids/ 1 
call  chlenldagdir, lendir) 
write(cform,90)  label, itau 
90  format(a8, i6.6) 
c 

filedg=  dagdi r(1: lendir)  //  cform 
c 

readlcdtg, • ( f 2 . 0 , f 6 . 0 ) • )  yr,  xdtg 
c 

c  open(uni t=44,f i le=f i ledg, form=' unformat ted' ) 
open(unit=44,f i le=f i ledg,form=‘formatted' ) 
c 

rd1=  idl 
rd2=  id2 
rlen=  len 

143  format(a64) 

write(44,144)  rdl , rd2, r len, yr, xdtg 
write(44,143)  capsun 
write(44,143)  (cap( i ) , i=1 , id2) 
write(44,144)  (x( i ) , i=1 , len) 

144  format(5e16.9) 
close(unit=44) 
print*,'  writing  1 , filedg 
print*,  capsun 

do  200  k=1 , id2 
print*,k,cap(k) 

200  continue 
c 

return 
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end 

subroutine  f ixitla, len, alow) 
c 

c  this  routine  is  used  to  adjust  certain  variables  which  during 
c  interpolation  are  given  unrealistic  values,  e.g.,  negative  vapor 
c  pressures, 
c 

c  *****************************************************************1 

c 

c 

implicit  real  (a-h,o-z) 
c 

real  a(len) 
c 

if (alow. ge. 0.0)  then 
do  10  i=1 , len 
a ( i ) =  max(a( i ) ,alow) 

10  continue 
else 

tem=  -alow 
do  20  i=1 , len 
a(i)=  minlal i ) , tem) 

20  continue 
endif 
return 
end 

subroutine  geotrpl im, jm, Im, lpout,psout,tt,phi ,sgeo,slp 
*,  plog, temps, phi ps) 
c 

c  this  routine  reads  in  the  model  forecasts  of  geopotential  and 
c  processes  them  so  they  can  be  output  as  standard  height  fields, 
c 

c  it  also  processes  the  temperature  output  fields  and  outputs 
c  these  as  standard  fields 
c 

dimension  plog(im*jm, lm+1 ),phips( im* jm, Ipout) 

*,  temps(im*jm, lpout),phi(im*jm. In) , tt( im* jm, lm+1 ) 

*,  phibot<im*jm),phisuflim*jm),sgeo(im*jm),slplim*jm) 

*,  psoutl Ipout) 
c 

c  convert  geopotentials  to  d-vatues  at  sigma  pressures 
c 

ograv=  1.0/9.80616 
do  11  k=1 , Im 
do  11  i=1,im*jm 
phi(i  ,k)=  phi(i,k)*ograv 
tt(i,k)=  tt( i ,k)-273. 15 

11  continue 
c 

p0log=  log(1000.0) 
do  13  i=1,im*jm 
phibot(i)=  0.0 
slplog=  log(slpli)) 

plog( i , lm+1 )=  maxis lplog,p0 log )+0. 0001 
13  continue 
c 

cal l  s2ptrp( im, jm, Im, Ipout, phi , phi bot.pl og, phi ps,psout , 0.5) 
c 

cal l  s2ptrp( im, jm, Im, Ipout , tt , tt( 1 , lm+1 ), plog, temps, psout, 0.5) 
c 

return 

end 

subroutine  glof lp( imax, jmax,glob,gf l ip) 
dimension  gt obi imax, jmax),gf l ipl jmax, imax) 
c 

do  10  j=1 , jmax 
do  10  i = 1 , imax 
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gflip(j,i)=  glob(i, jmax+1-j) 

10  continue 
return 
end 

subroutine  hadly(im, jm, tm, jtrun, mlmax, weight, poly 
*,  msort, stream) 

dimension  meq0( jtrun), msort (mlmax), spcoef (mlmax, 2) 

*,  poly(mlmax*2, jm) .weigh t( jm) , work4( j trun) 

*,  work5(im, jm),stream(jm,lm) 

m=  0 

do  150  ml=1 , mlmax 
if(msort(ml).ne.1)  go  to  150 
m=  m»1 
meq0(m)=  ml 

work4(m)=  m/sqrt(4.0*m*m-1 .0) 

150  continue 

do  450  k=1 , Im 

do  452  j=1 , jm 
do  452  i=1 ,  im 
work5(i,j)=  stream(j.k) 

452  continue 

call  tranrs(0,1 , poly, weight, works, spcoef) 
work5( 1 , 1 )=0.0 

work5(2, 1 )=0.5*spcoef (meq0(3), 1 )*work4(2) 
jtrun2=jtrun-2 

do  446  mm=2,jtrun2 
ml=  meqO(inn) 
ml  2=  meq0(nm*2) 

work5(mm*1 , 1 )=(nm*spcoef (ml  2, 1 )*work4(mm+1 )- (mm*1 )*work4(mm) 
*  *spcoef  (ml ,  1  ))/(rrm*(rmn-1 ) ) 

446  continue 

ml1=  meq0( jtrun-1) 

workS ( jtrun, 1 )=- spcoef (ml  1 , 1 )*work4( j  trun)/( j trun- 1.0) 
call  zi lch(spcoef ,mlmax*2) 

do  447  1=2, jtrun 

spcoef (meq0( l ), 1 )=  workS (1,1) 

447  continue 

cat l  transr(1 , poly, spcoef ,work5, 0,0,0) 

do  460  j=1 , jm 
stream(j,k)=  work5(1,j) 

460  continue 

450  continue 


return 

end 

subroutine  hterp( ik, z, im, jm,glob,xin,yin, ix, jy.alat , v) 

dimension  glob( ix, jy), wl ( im, jm+2) 

*,  alat( jm*2),xin(ix, jy),yin(ix, jy) 

real  v(1000) 


c 

dimension  z( im, jm) 
c 

save  lenx 
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c 

c 


data  rad/6. 375e6/ 


pi  =  4.0*atan(1 .0) 
pio2=  pi*0.5 
if(ik.eq.O)  then 
lenx=  ix*jy 
tem=  im/f loat(ix) 
do  10  i=2,ix 

xin(i,1)=  1.0+tem*(i-1.0) 

10  continue 

xinfl, 1 )=  1.00001 
do  15  j=2, jy 
do  15  i=1,ix 
15  xin(i,j)=  xin(i,1) 
c 

jym=  jy- 1 
do  17  j=2, jym 
tem=  -pio2+(j-1 .0)*pi/jym 
do  17  i=1,ix 
17  yin(i,j)=  tem 
do  12  1=1 , ix 
yin(i,1)=  -pio2 
yin(i.jy)=  pio2 
12  continue 
endif 
c 

do  22  i=1,im*jm 
22  w1<i,2)=  z(i,1) 
call  sumit(im,z,ps) 
call  sumit(im,z(1, jm),pn) 
do  25  i=1 , im 
w1(i,1)=  ps/im 
25  w1(i,jim-2)=  pn/im 

cal  l  bcubiy(ik,2,w1 ,  im,  jm*-2,glob,  lenx.alat ,xin,yin,v,  1 .0, 1 .0) 
c 

return 

end 

subroutine  lablslkid, Irec, iplev, isize, idtg, idtgs, idtge.gmean 
*,  capsun, ci ,co,add,amp, l lab, ident) 
c 

character*192  ident 
character*8  idtg, idtgs, idtge 
character*24  capsun 
character*^  iplev.ctau 
character*1  if  lap 
character*3  Irec 
c 

character*24  (lab 
c 

k1=  1 

do  12  k=1 ,24 
ident<k1:k1*7)=' 
k1=  kl+8 
12  continue 
c 

ident(169:192)=  capsun 
if(iplev.eq.'  ')  then 
ident(49:72)=  l lab 
ident(73:88)=  '  ' 

ctau=  '  O' 
else 

wri te( ident (49: 88), 100)  l lab, iplev 
100  format(a24, '  AT  \a4,'  MBS') 
ctau=  iplev 
endif 
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write(ident(89:136),200)  idtgs, idtge.gmean 
200  formate 'TIHEAVE  FROH  1 ,A8, 1  TO  ',A8, '  mean=',f9.3) 
c 

iftap=  '$' 
lenx=  10536 

if(isize.gt.O)  lenx=  24+288*145 
if(isize.gt.O)  iflap='N* 

ident(1 :16)  =  lrec  //  if  lap  //  ctau  //  idtgs 
ident(48:48)='0' 
write(ident(17:24),210)  lenx 
210  format(2x, i6) 

ident(33:40)='  T I HEAVE' 
c 

if (ci . It. 0.01 )  then 

write(ident(137:168>,  '(4f8.4)')  ci , co, add, amp 
c 

else  if (abs(add).gt. 100.0)  then 
uri te( i dent (137: 168), ' <4 f 8 . 1 ) ' )  ci ,co, add, amp 
c 

else 

uri  ted  dent  (137: 168), 1  (4f8.3) ' )  ci  ,co,add,amp 
c 

endi  f 
c 

print  300,  ident 
300  format(24a8) 
return 
end 

subroutine  levget( ip, pout, Ipout.klev) 
dimension  pout (l pout) 
c 

klev=  0 

do  10  k=1,lpout 

i f ( ip.eq. int(pout(k) ) )  then 

klev=  k 

return 

endi  f 

10  continue 
return 
end 

subroutine  p2kapa(imjm, lm,asig,bsig,pt,psigc,plog,pk2,pl l ) 
c 

dimension  asig( lm+1 ),bsig( lm+1 ),pt( imjm),pk2( imjm, Im) 

*,  psigc( imjm, lm),plog( imjm, lm),pl l ( imjm, tm) 
c 

c  compute  pressures  on  the  sigma  surfaces  of  the  forecast  model 
c  and  raise  to  the  capa  power, 
c 

dcapa=  3.5 
capa=  1.0/3. 5 
capap1=  1.0+capa 

p0=  1000.0 

opok=  1.0/p0**capa 
p0log=  log(pO) 
tem1=  opok 
c 

px=  1.0/capapl 
do  19  i=  1,imjm 
pwrk2=  asig(2)+bsig(2)*pt( i ) 
pk2(i,1)=  opok*pwrk2**capa 

psigc( i , 1 )=  (pk2( i , 1 )*pwrk2-tem1 )*px/(pwrk2- 1 .0) 
plog( i ,1 )=  p0log+dcapa*log(psigc( i , 1 )) 
pll(i,1)=  exp(plog( i , 1 )) 

19  continue 
c 

do  28  k=2, Im 
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do  27  i*1 , imjm 
pwrk1=  asig(k)*bsig(k)*pt( i ) 
pwrk2=  asig(k*1)+bsig(k+1)*pt(i) 
pk2(i,k)=  opok*pwrk2**capa 

psigc(i,k)=(pk2(i,k)*p«rk2-pk2(i ,k-1)*p«rk1)*px 

*  /(pwrk2-pwrk1 ) 

plog( i ,k)=  pOlog+dcapa*log(psigc( i ,k)) 
pll(i,k)=  exp(plog(i,k)) 

27  continue 
2ft  continue 

return 

end 

subroutine  ritasc(ihdg, len,x) 
dimension  x(len) 
logical  first 

character*192  ihdg 
data  first/. true./ 

if(first)  then 
f i rst=. false. 

openluni t=16, f i le='/scr/rosmond/pcmdi/pcmf tds/asf i  te1 

*  ,forro=' formatted' ) 
endif 

write(16,100)  len.ihdg 
write(16,200)  (x( i ), i=1 , len) 

100  format( i8,a192) 

200  format(8e10.4) 
return 
end 

subroutine  s2pcs( im, jm, lm,y,boty,pin,dout,pout,fac) 

dimension  y( im, jm, lm*1),pout( lm),boty(im, jm),dout( jm, Im) 
*,  pin(im,  jm,  lim-1) 

dimension  yh( im, lm+1 ),pinh( im, lm*1),doh( im, Im) 

*,  poh(im, Im) 

eps=  1.0e-6 
ten=  0.5 
imlm=  im*lm 
do  5  k=1,lm 
do  5  i=1,im 
poh( i ,k)=  pout ( k ) 

5  continue 

do  1  j=1,jm 
do  4  i=1 , im 
yh(i,lm*1)=  boty(i,j) 

4  continue 
do  2  k=1 , Im 
do  2  i=1 , im 
yh( i ,k)=  y(i,j,k) 

2  continue 

do  3  k=1 , lm*1 
do  3  i=1, im 

3  pinh( i ,k)=  pin(i,j,k) 

cal l  bcubv(0,yh, im, lm»1 ,doh, imlm.pinh.poh, ten) 
xfac=  fac/im 
do  6  k*1,lm 
dout(j.k)*  0.0 
do  6  i=1,im 

<fout(j,k)=  dout(  j  ,k)*doh(  i  ,k)*xfac 

6  continue 
1  continue 

return 
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c 


c 


c 


c 


c 

c 

c 

c 

c 


c 


end 

subroutine  s2ptrp( im, jm, Im, lpout,y,boty,pin,dout, pout, ten) 

dimension  y(im, jm, lm+1 ), pout ( Ipout ),boty(im, jm),dout( im, jm, l pout ) 
*,  pin(im, jm, lm+1) 

dimension  yh(im,  lm*1  ),pinh(  im,  lm*-1  ),doh(  im,  l pout) 

*,  poh(im, Ipout) 

eps=  1.0e-6 
imlm=  im*lpout 
do  5  k=1, Ipout 
do  5  i=1,im 
poh(i ,k)=  pout(k) 

5  continue 

do  1  j=1,jm 
do  4  i*1,im 
yh( i , lm»1 )=  boty(i,j) 

4  continue 
do  2  k=1,  Im 
do  2  i=1, im 
yh(i,k)=  y(i, j,k) 

2  continue 

do  3  k=1 ,  Iith-1 
do  3  i=1,im 

3  pinh(i,k)=  pin(i,j,k) 

cal  l  bcubv(0,yti,  im,  lnn-1  ,dch,  imlm,pinh,poh,  ten) 

do  6  k=1 , Ipout 

do  6  i=1 , im 

dout( i , j ,k)=  doh( i ,k) 

6  continue 
1  continue 

return 

end 

subroutine  sdet 

parameter  (im=144,  jm=im/2,  lm=  18,  lpout=6) 

parameter  (jtrun=  2*( (1+( im-1 )/3)/2) ,  mlmax=  < j trun+1 )* j trun/2) 

parameter  (idim2=  2*mlmax) 

parameter  (nspec=4*lm+1 ) 

parameter  (numave=  40+21*lpout+2+14*lpout+7) 
parameter  (imjm=im*jm) 


parameter  (lccn=  350+lm*(12+5*lm)) 


dimension  ccn(lccn) 
equivalence  (c,ccn) 
character's  idtg 

common/cnstnt/  c(300), itau, ithist, idtg(3), sarray( Im, Im) 
*,  carray(lm, lm),evecin(lm, lm),evectr(lm, lm),tmcor(lm, Im) 
*,  eigvaK lm),spalm( lm),pmcor( lm),pbcor( Im), tmeanC Im) 

*,  asigl lm+1 ),bsig( lm+1 ),dasig( lm),dbsig( Im) 

*,  sige( lm*1 ),dsig( lm),pkfac( Im) 

*,  ptmean,ccpad(100) 


equivalence 
1  (c(1),jmm1), 

2,  (c(5),iday), 

3,  <c(7),itauo), 

4,  (c(10), jsplit), 

5,  (c(13),taue), 


(c(3), jmm2), 
(c(6),tofday) 
(c(8), julian), 
<  c ( 1 1 ) , taui ), 
(c(14),tice). 


(c(4), Increc) 

<c(9), jskip) 

(c(12),dtrad) 

(c(15),forwrd) 
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c 


c 


c 

c 

c 

c 

c 

c 


c 


6,  <c(16),cosd), 

7,  (c(19),eccld), 

8,  (c(22),prevap), 

9,  <c<25), omega), 

1,  (c(28),pi2), 

2,  (c(31),cp), 

3,  (e(41),f luxcl), 

4,  (c(44),tauh), 

5,  <c(47),mombdg), 

6,  (c(50),p00), 

7,  (c(53),hybrd), 

8,  (c(56),dgtime), 

equivalence 
1  (c(68),taud), 

2,  (c(70),pcon), 

3,  <c(73),skew), 

4,  (c(76), ithfg), 

5,  (c(79),tseres), 

6,  (c(93),daypyr), 

7,  (c(96), iqnx), 

8,  (c(99),eccn), 

9,  (c(102),critl), 

1,  (c(105),pcp), 

2,  (c(109), ipfcst), 

3,  (c(112), ipoutp), 

4,  (c(115), icnmax), 

5,  (c(118),nmodev), 

6,  (c(121),uvmax), 
8,  <c(135l,nrnt). 


<c(17),ddraf t), 

<c(20),rbear), 

(c(23),radsq) 

(c(26), tauave), 

<c(29),rgas), 

(c(32),hltm), 

(c(42), isunfx), 

(c(45),krnit), 

(c(48),evaprh), 

(c(51),dta), 

(c(54),ksgeo), 

(c(57),ecadv) 


(c(69),dt) 

(c(71),pconkt), 

<c(74),grav2), 

(c(77), itaufg), 

(c(80),prftop), 

(c(94),day), 

(c(97),perhl), 

(c(100),rad), 

(c(103),p608) , 

(c(106),frqdib), 

(cCIIO), ipsigp), 

(c(113),tau), 

(c(116),frqlim), 

(c(119), dinit), 

(c(123),facrad), 

(c( 136) , capapl ) , 


(c< 18) , jtrunp) 
(c(21),prnti j) 

(c(27),pi ) 
(c(30),capa) 
(c(33),fim) 
(c(43) ,ptop) 
(c(46), Is) 
(c(49), vis) 
(c(52),ops) 
(c(55),pok) 


(c(72) , tsl ) 
(c(75),delmf ) 
(c(78), Isx) 
(c(81),dtgpcm) 
(c(95),rotper) 
<c(98),decmax) 

<  c< 1 01 ) ,grav) 
(c(104), Imid) 
(c(107),cdtgoi ) 
(c(111), ipinit) 
(c(114),drgdry) 
(c(117),beta) 
Cc( 120) , taudb) 
(c(134),gamfac) 
(c(137),cuprh) 


equivalence 
1  (c(138),novar), 

2,  (c(  141), nozone), 

3,  (c(144), incrup), 

4,  (c<151),nsdedy), 

5, 


(c(139),cdet), 
(c(142),sstclm), 
Cc(145), incwnd), 
(c(152),kmonth), 
(c(157),rsdist), 


(c(140),dcapa) 

(c(143),updayc) 

(c(146),pcmdi) 

(c(153),mnthdy) 

(c(158),sind) 


******************************************************************** 


logical  first 

save  sob,perhlr,reccn,h,maxday,y, i iday 
data  first/. true./ 

f(x)  =  reccn*(1 .-eccn*cos(x-perhlr))**2 

if(first)  then 
iqnx=  80 

pib180  =  pi/180. 

sob  =  sin(decmax*pib180) 

perhlr  =  perhl  *  pib180 

reccn  =  1 ./(I . -eccn*eccn)**1 .5 

h  =  pi 2/daypyr 

maxday  =  daypyr  ♦  .1 

endif 


c 

c  set  the  date 

c 

call  daynuiK idtg(2), idy,kmonth,nnthdy,mh,nsdedy,mhy) 
c 

c  integrate  for  the  earth's  position 

c 


if  (.not. first  .and.  nsdedy  ,ne.  iqnx)  go  to  105 
first  =  .false, 
y  «  0. 
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* 


i i day  =  nsdedy  -  i qnx 
if  (iiday  .eq.  0)  go  to  110 
if  (iiday  .It.  0)  iiday  =  iiday  +  maxday 
100  iiday  =  iiday  -  1 
105  tl  =  f(y)*h 

t 2  =  f (y+t1*.5)*h 
t3  =  f (y+t2*.5)*h 
t4  =  f(y+t3)*h 

y  =  y  +  (t1+2.*(t2+t3)+t4)/6. 
if  (iiday  .gt.  0)  go  to  100 
c 

c  compute  declination  and  earth-sun  distance 
c 

110  sind  =  sob*sin(y) 

cosd  =  sqrt(1 .-sind*sind) 

rsdist  =  ((1.-eccn*eccn)/(1.-eccn*cos(y-perhlr)))**2 
c 

print  120,  nsdedy, kmonth.mnthdy, sind 
120  formate  nsdedy=' ,  i3, ' ,  month=* ,  i  2, 1 ,  day=‘,i2 
sine  of  solar  declination=' ,f6.4) 


c 


c 

c 

c 

c 

c 


c 


c 

c 

c 


return 

end 

subroutine  solar(s0cosz, idtgx.sint .cost) 
parameter  (im=144,  jm=im/2,  tm=  18,  lpout=6) 

parameter  (jtrun=  2*((1+( im- 1 )/3)/2),  mlmax=  ( j trun+1 )* j trun/2) 
parameter  (idim2=  2*mlmax) 
parameter  (nspec=4*lm+1 ) 

parameter  (nunave=  40+21*lpout+2+14*lpout+7) 
parameter  (imjm=im*jm) 


parameter  (lccn=  350+lm*(12+5*lm)) 

******************************************************************** 


dimension  ccn(lccn) 
equivalence  (c,ccn) 
character*8  idtg 

conmon/cnstnt/  c(300), itau, i thist , idtg(3) ,sarray( Im, Im) 

*,  carray( Im, lm),evecin(lm, lm),evectr( Im, Im), tmcor( Im, Im) 

*,  eigval( lm),spalm( lm),pmcor( lm),pbcor( Im), tmean( Im) 

*,  asig(lm+1),bsig(lm+1),dasig(lm),dbsig(lm) 

*,  sige(lm+1),dsig(lm),pkfac(lm) 

*,  ptmean,ccpad(100) 

******************************************************************* 


equivalence 
1  (c(1),jmm1), 

2,  (c(5),iday), 

3,  (c(7), i tauo), 

4,  (c(10), jsptit), 

5,  (c(13),taue), 

6,  (c(16),cosd), 

7,  (c(19),eccld), 

8,  (c(22),prevap), 

9,  (c(25), omega), 

1,  (c(28),pi2), 

2,  (c(31 ),cp), 

3,  (c(41),f luxcl), 

4,  (c(44),tauh), 

5,  (c(47),mombdg), 

6,  (c(50),p00), 

7,  (c(53),hybrd), 


(c(3),  jrrm2), 

(c(6),tofday) 

(c(8>, julian), 

(c(11),taui), 

(c(14),tice), 

(c(17),ddraft), 

(c(20),rbear), 

(c(23),radsq) 

(c(26),tauave), 

(c(29),rgas), 

(c(32),hltm), 

(c(42), isunfx), 

(c(45),krnit), 

(c(48),evaprh), 

(c(51),dta), 

(c(54),ksgeo). 


(c(4),tncrec) 

(c(9), jskip) 
(c(12),dtrad) 
(c(15), forwrd) 
(c( 18) , jtrunp) 
(c(21 ),prnti j ) 


(c(27),pi ) 
(c(30),capa) 
(c(33),f im) 
(c(43),pt op) 
(c(46) , Is) 
(c(49),vis) 
(c(52),ops) 
(c(55),pok) 
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c 


c 


c 

c 

c 

c 


c 

c 

c 

c 

c 


8,  (c(56),dgtime), 

equivalence 
1  (c(68) , taud) , 

2,  (c(70),pcon), 

3,  <c(73),skew), 

4,  <c(76), ithfg), 

5,  <c(79),tseres), 

6,  (c(93),daypyr), 

7,  (c(96),iqnx), 

8,  (c(99),eccn), 

9,  <c(102),critl), 

1,  (c(105),pcp), 

2,  <c(109),ipfcst), 

3,  (c(112), ipoutp), 

4,  <c(115), icnmax), 

5,  (c(118),nmodev), 

6,  (c(121),uvmax), 
8,  (c(135),prnt), 

equivalence 
1  (c(138),novar), 

2,  (c( 141), nozone), 

3,  (c(144), incrup), 

4,  (c(151),nsdedy), 

5, 


(c(57),ecadv) 


(c(69),dt) 

(c(71),pconkt), 

(c(74),grav2), 

(c(77), itaufg), 

(c<80),prf top), 

<c(94),day), 

(c(97),perhl), 

<c(100),rad), 

(c(103),p608), 

(c(106),frqdib), 

(c(110), ipsigp), 

(c(113), tau), 

(c(116),frqlim), 

(c(119), dinit), 

(c<123),facrad), 

(c(136),capap1), 


(c(139),cdet), 
(c(142),sstclm), 
(c(145), incwnd), 
(c(152),kmonth), 
(c(157),rsdist). 


(c(72), Isl) 

(c(75),delmf ) 

(c(78), Isx) 

(c(81),dtgpcm) 

(c(95) , rotper) 

(c(98),decmax) 

(c(101),grav) 

(c(104), Imid) 

(c(107),cdtgoi ) 

<c( 111 ), ipini t) 

(c(114),drgdry) 

(c(117),beta) 

(c(120),taudb) 

(c(134),gamfac) 

(c(137),cuprh) 


(c(140),dcapa) 

(c(143),updayc) 

(c(146),pcmdi) 

(c(153),mnthdy) 

(c(158),sind) 


************* 


character*8  idtgx 

dimension  sOcosz(imjm),sinl( jm),cosl( jm) 
data  solcst/1365.0/,  pi/3.1415926/ 
idtg(2)=  idtgx 

compute  siderial  day  and  solar  declination  for  this  day 
call  sdet 


c  adjust  solar  constant  for  this  day;  see  Paltridge  and  Platt  (1976) 
c  use  value  at  1200  for  daily  mean 
c 

fday=  nsdedy-0.5 

tho=  min(2.0*pi ,2.0*pi*fday/365.0) 

tem2=  2.0*tho 

sin0=  sin(tho) 

cos0=  cos(tho) 

sin20=  sin(tem2) 

cos20=  cos(tem2) 

s0=  so lest* ( 1 . 0001 1 0+0 . 034221 *cos0+0 . 00 1 280*s i nO 
*  +0.000719*cos20+0.000077*sin20) 
c 

c  integral  of  cos(Z)  over  daylight  period,  multiplied  by  sO 
c 

do  10  i=1,imjm 
j=  1+(i-1)/im 
a=  sinl(j)*sind 
b=  cosl(j)*cosd 
arg=  a/b 

if(abs(arg).le.1.0)  then 
h0=  acos(-arg) 
s0cosz(i)=  a+b/h0*sin(h0) 
s0cosz(i)=  s0cosz(i)*h0/pi 
sOcosz(i)-  s0*s0cosz(i) 
else  if(arg.gt.l.O)  then 
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s0cosz(i)=  sO*a 
else  if(arg.lt.-I.O)  then 
sOcosz(i)=  0.0 
end  if 

c  if(mod(i,144).eq.O)  print  100,  j ,a,b,arg,h0,s0cosz( i ) 
c  100  formatOx, '  j  =  ' , i3, '  a=',e10.4,'  b=',e10.4,1  arg=',e10.4 
c  hC=',e10.4,‘  sOcosz=‘ ,e10.4) 

10  continue 
return 
end 

subroutine  suncos( im, jm.cosl ,x,gmean) 
c 

dimension  cosl ( jm),x( im, jm) 
logical  first 
c 

save  denom 
c 

data  first/. true./ 
if (first)  then 
dertom=  0.0 
do  1  j=1,jm 

1  denom=  denom+cos l ( j ) 
denom=  1 .0/( im*denom) 
first=. false. 

endif 

c 

gmean=  0.0 
do  2  j=1 , jm 
do  2  i=1,im 

2  gmean=  gmean+cos l ( j )*x( i , j ) 
gmean=  gmean*denom 

return 

end 

subroutine  svar( imjm, Im, kmid.pt ,sgeo,pdi f f , tsave.phi , tt, tmppl ,slp) 
c 

c  subroutine  to  process  single  level  surface  variable  output 
c 

c  this  routine  reads  the  model  fields  of  pbl  variables  and  ground 
c  hydrology  variables  and  processes  them  so  they  can  be  output  as 
c  standard  fields, 
c 

dimension  sgeo( imjm) ,pt( imjm) ,pdiff( imjm), tsave(imjm) 

*,  tt( imjm, lm*1 ), tmppl (imjm), phi (imjm, lm),slp( imjm) 
dimension  ppp(imjm) 
c 

c  delta  p  method  for  finding  sea  level  pressure 
c 

c  add  pdiff  to  terrain  pressure  to  get  forecast  sea  level  pressure 
c  compute  correction  to  pdiff  due  to  change  in  temperature  in 
c  bottom  half  of  atmosphere  during  forecast 
c 

cp=  1005.45 
rgas=  cp/3.5 

temx=  1 .0/(rgas*log(10.0)) 
ocp=  1.0/cp 
c 

c  compute  tIOOO  from  deep  layer  thickness  change  and  adiabatically 
c  from  bottom  model  level  pressure. 

c  we  constrain  the  1000  mb  temperature  to  be  between  isothermal 
c  and  adiabatic  w.r.t.  the  bottom  model  level  temperature 
c 

do  40  i=1 , imjm 

ppp( i )=  temx*(phi ( i ,kmid)-sgeo( i ))-tsave( i ) 

p2=  tmppl (i)+ppp(i) 

p1=  tt(i,lm)+ocp*sgeo(i) 

p2=  min(p2,p1) 

p2=  max(p2,tt(i,lm)) 
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c 

c 

c 


c 


c 


c 


c 

c 


c 

c 


adjust  pdiff  due  to  change  in  subterrain  temperature 

p1  =  p2-tmpp1(i) 
tt(i,lm+1)=  p2 

p2=  rgas*tt(i,lm«-1)*tt(i,lnHl) 
p2=  sgeo( i >*(pt( i )+pdi f f ( i ) )/p2 
ppp(i)=  p1*p2 

slp(i)=  pt(i)*pdiff(i)-p1*p2 
40  continue 

return 

end 

subroutine  trngra( l l , cim, poly, dpoly,s, dipt ,dtpl ) 

parameter  (im=144,  lm=18,  lpout=6) 
parameter  (jm=  im/2,  im2=  im/2,  jm2=  jm/2) 

parameter  (jtrun=  2*((1+( im-1 )/3)/2) ,  mlmax=  ( j trun+1 )* j trun/2) 
parameter  (imlm=im*lm,  imjm=im*jm,  idim2=  mlmax*2) 

parameter  (mlx=  ( jtrun/2)*(( jtrun+1 )/2) ) 
parameter  (jump=  iim-3) 

common/  fftcom/  trigsC im,2), ifax(19) 

dimension  cim(mlmax) ,poly( idim2, jm2),dpoly( idim2,  jm2) 

*,  s(mlmax,2, lm),dlpl( im, jm, lm),dtpl(im, jm, Im) 

dimension  cu( jump, jm),cv( jump, jm),work( imjm,2) 

do  15  k=1,U 
do  27  m=  1,jump*jm2 
cu{m,1)=  0.0 
cu(m, jm2+1 )=  0.0 
cv(m, 1)=  0.0 
cv(m, jm2+1)=  0.0 
27  continue 


do  20  j=1,jm2 
jj=  jm+1-j 
c 

ml=  2*mlx 
CDIR3  IVDEP 


do  50  m=2, jtrun 
ml=  ml+1 
mm=  2*m-1 
mp=  um*1 
cu(mm, j )= 

2 

+cim(ml )*poly(ml , j )*s(mi ,2,k) 

cu(mm, j  j )= 

+cim(ml )*poly(mt , j )*s(ml ,2,k) 

cu(mp, j )= 

-cim(ml)*poly(ml , j }*s(ml , 1 ,k) 

cu(mp, j  j  >= 

-cim(ml )*poly(ml , j )*s(mi , 1 ,k) 

cv(mn,  j  )= 

-dpoly(ml , j )*s(ml , 1 , k) 

cvCmm, jj)= 

♦dpoly(ml, j)*s(ml,1,k) 

cv(mp, j )= 

-dpoly(ml, j)*s<ml,2,k) 

cv(mp, jj)= 

♦dpoly(ml, j)*s(mt,2,k) 

50  continue 

m1=  0 

do  30  l=jtrun-1,1,-2 
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« 


4 


COIRS  1VDEP 

do  40  m=1 , l 
mm=  2*m-1 
mp=  rmw-1 
ml=  m*m1 
mk=  ml+mlx 


c 

cu(mm,  j  )=  cu(mm,  j  )+cim(ml  )*poly(ml ,  j  )*s(ml ,2,k) 

*  +cim(mk)*poly(mk, j)*s(mk,2,k) 
c 

cu(mm,  j  j  )=  cu(mm,  j  j  )+cim(ml  )*poly(ml , j )*s(ml , 2, k) 

*  -cim(mk)*poly(mk, j)*s(mk,2,k) 
c 

cu(mp,  j)=  cu(mp, j)-cim(ml)*poly(ml, j)*s(ml,1,k) 

*  -cim(mk)*poly(mk, j)*s(mk,1,k) 
c 

cu(mp, j  j )=  cu(mp, j  j )-cim(ml )*poly(ml , j )*s(ml , 1 , k) 

*  ♦cim(mk)*poly(mk, j )*s(mk, 1 ,k) 
c 

cv(mm,  j  )=  cv(mm,  j )-dpoly(ml ,  j  )*s(ml ,  1 ,  k) 

*  -dpoly(mk, j)*s(mk,1,k) 
c 

cv(mm, j  j )=  cv(mm, j  j )+dpoly(ml , j )*s(ml , 1 ,k) 

*  -dpolyfmk, j)*s(mk,1,k) 
c 

cv(mp, j )=  cv(mp, j )-dpoly(ml , j )*s<ml ,2,k) 

*  -dpolyCmk, j)*s(mk,2,k) 
c 

cv(mp, j  j )=  cv(mp, j  j )+dpoly(ml , j )*s(ml ,2, k) 

*  -dpoly(mk, j)*s(mk,2,k) 
c 

40  continue 


m1=  ml  +  l 
30  continue 
20  continue 
c 

call  rf ftmlt(cu, work, tri gs, i fax, 1, jump, im, jm, 1) 
cal l  rfftmlt(cv,work, trigs, ifax, 1 .jump, im, jm, 1 ) 
c 

do  22  j=1,jm 
CDIR3  IVDEP 

do  22  i =1 , im 
dlpUi,  j,k)=  -cu(i,j) 
dtpl(i,j,k)=  -cv(i, j) 

22  continue 
c 

15  continue 
return 
end 

subroutine  var(imjm,fac,x,xv,avsum) 
c 

dimension  x( imjm),xv( imjm) ,avsum( imjm) 
c 

do  10  i=1,imjm 

avsum( i )  =  fac*xv( i )-(x( i )*fac)**2 
10  continue 
return 
end 

subroutine  wndf ix(x, imx, jinx) 
c 

c  programmer  --  rosmond  t.  (02  nov  1989) 
c 

c  this  routine  computes  values  of  wind  field  at  poles  and 
c  de-weights  winds  with  cos(lat) 
c 

c  ******************************************************************** 
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c 

dimension  x(imx,jmx) 
dimension  cosl(200) 
save  cost, jmm,im2,pi,jmsav 
c 

data  jmsav/-100/ 
c 

ifCjmx.ne. jmsav)  then 
pi=  4.0*atan<1 .0) 
do  3  j=1,jmx 

cost( j)=  pi*(0.5-( j-1.0)/< jmx-1.0)) 
cost(j)=  cos(cosl(j)) 
cosl(j)=  6.375e6/cost( j ) 

3  continue 
jmnt=  jmx-1 
im2=  imx/2 
jmsav=  jmx 
endif 
c 

do  1  j  =2 ,  jmro 
do  1  i=1, imx 
x(i,j)=  x(i,j)*cosl(j) 

1  continue 
c 

CD  IRS  IVDEP 

do  5  i=1 , im2 

x(i,1)=  0.5*(x(i,2)-x(i+im2,2)) 
x(i+im2,1)=  -x(i,1) 
x(i, jmx)=  0.5*(x(i, jmm)-x( i+im2, jmm)) 
x(i+im2, jmx)=  -x(i,jmx) 

5  continue 
return 
end 

subrout i ne  womegaC im, jm,  Im, onocos ,  b sig, das i g, dbs i g,  pt 
*,  dtpl,dlpl,rdiv,ut,vt,sd) 
c 

dimension  bsig( lm*1 ),ut( im* jm, Im) , vt ( im* jm, Im) ,onocos( im* jm) 
*,  dlpl( im*jm),dtpl(im*jm),pten( im*jm),rdi v( im*jm, Im) 

*,  dasigt lm),dbsig( lm),pt( im*jm) 
c 

dimension  g( im* jm, Im) ,sd( im* jm, Im) 
c 

k=  1 

do  22  i=1,im*jm 

g(i,k)=ut(i,k)*dlpl(i)*onocos(i)+vt(i,k)*dtpl(i) 
pten( i )=g( i ,k)+rdi v( i ,k)*pt( i ) 
pten( i )=  -pten(i )*dbsig(k)-rdiv(i ,k)*dasig(k) 
sd( i ,k+1 )=  pten( i ) 
sd( i , 1 )=  0.0 
22  continue 
c 

do  2  k=2, lm-1 
do  2  i=1,im*jm 

g( i ,k)=  ut( i ,k)*dlpl( i )*onocos( i )+vt( i ,k)*dtpt ( i ) 

8=  g(i,k)+rdiv(i,k)*pt(i) 

pten< i )=  pten(i )-a*dbsig(k)-rdiv< i ,k)*dasig(k) 
sd(i,k+1)=  pten(i) 

2  continue 
c 

k=  Im 

do  24  i=1,im*jm 

g( i ,k)=  ut C i ,k)*dtpl ( i )*onocos( i )+vt( i , k)*dtpl ( i ) 
a=  g(i ,k)+rdi v( i ,k)*pt( i ) 

pten(i )=  pten( i )-a*dbsig(k)-rdi v( i ,k)*dasig(k) 

24  continue 
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c  vertical  velocity 
c 

do  3  k=2, Im 
do  3  i=1,im*jm 

sd(i,k)=  sd( i ,k)-bsig(k)*pten( i ) 

3  continue 
c 

do  450  k=1 , Im 
c 

if(k.lt.lni)  then 
c 

do  443  i=1,im*jm 

sd( i ,k)=bsig(k+1 )*(0.5*(g( i , k)+g< i , k+1 ) )+pten( i ))+sd( i ,k+1 ) 
443  cont i nue 

c 

else 

c 

do  441  i=1,im*jm 

sd(i,lm)=  g(i,lm)+pten(i) 

441  continue 
end  if 

450  continue 
c 

return 

end 
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